NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


> 

THESIS  A 

OPTIMIZATION  OF  AN  INTERNALLY 
FINNED  ROTATING  HEAT  PIPE, 

by  '/  j 

i  / 

I 

William  A.  Davis,  Jr[ 

j  Septembei*-“19  8  0 

/  . 

Thesis  Advisor:  P.  J.  Marto 

Approved  for  public  release;  distribution  unlimited 


81  3  16  079 


$CCU*«TV  CLASSIFICATION  Of  THIS  • 


REPORT  DOCUMENTATION  PACE 


4.  TITLE  luFlIllal 


Optimization  of  an  Internally 
Finned  Rotating  Heat  Pipe 


RKAO  INSTRUCTIONS 

•BFORK  COMRLETOfC  FORM 


RCCIFtCNT'S  CATALOG  NUMSCR 


»■  TVA*  OF  »I»04T  »  P(MOO  COVlRtO 

Master's  Thesis; 
September  1980 

I.  FIRFORMINO  OHO.  RSFORT  numRCR 


7.  AuTmOF I  mi 


William  A.  Davis,  Jr. 


»  FCRFORMIN4  OROANI2ATION  M  AMI  AND  AOORCI1 

Naval  Postgraduate  School 
Monterey,  California  93940 


FRQORAM  (LEM (NT.  FFOjECT  T  »»k 
ARIA  *  WORK  UMtT  NUMIIRI 


>  »  COMTROLUNO  OfflCl  H  AM  E  *N0  A0OHCSS 

Naval  Postgraduate  School 
Monterey,  California  93940 


I 


MOnTTOHinG  AGENCY  name  A  AOO«C$ HU  t 


Cmntrmttlmg  Ottimm) 


I#.  OlftTftlftuTlON  ITATIMIMT  (ml  thlm  RaN#D 


12.  ftCPOIAT  OATS 

September  1980 


NUMRCR  OF  FADES 

106 


il.  ICCumitv  cl  All.  (mi  m it 

UNCLASSIFIED 


Approved  for  public  release;  distribution  unlimited 


17.  OlSTUi  •uTlOM  STATCMfMT  f ml  IA«  mrnmttmm  t  rnmtmtmxl  It i  fl#«*  JO,  Il  UllafWl  fmm  Mm*mrf) 


19  X  €  y  VOROI  <C «*  rmmmrmm  t**m  it  rtmmmmmmrf  mn4  Ifnitty  Wf«A  wAfO 

■  Heat  Pipe;  Rotating  Heat  Pipe;  Optimization 


10  A  It  TRACT  (Cmmilmtm  an  raaaraa  ,IM  if  naaaaaaar  an!  llaniifr  *T  Waal  — >«) 

A  finite  element  formulation  was  used  to  solve  the  steady-state 
two-dimensional  conduction  heat  transfer  equation  in  the  condenser 
wall  section  of  an  internally  finned  rotating  heat  pipe.  A  FORTRAN 
program  using  this  method  was  coupled  with  the  C0PES/C0NMIN  program 
for  sensitivity  analysis  of  design  variables  and  for  automated 
design  of  the  internal  heat  pipe  geometry. 

With  water  as  the  working  fluid,  numerical  results  obtained  for 
coD£er_^nd_^ai£lej^_i>^o®^>Jigi^<i^^_>£ond£n^£r_<sections_iindicateid__i 

DO  rjSTn  1473  coition  of  '  NOV  II  il  OMOLCTI 

(Pace  1)  i/N  iitoi*ii4- ami  —  i  ..I,.  - _ -  .  — y 

'  .  UtuFfv  Cl  AFAlFlC  ATiOn  OF  TniI  Fa«I  i'RNni  Df«  Iaiax 


•ccufity  cl  aiii  fic  at  ion  of  tnii  ra««  i 


Oata  RaiaraF) 


amiimo'igm.rii  nsasan 


that  for  the  maximum  heat  transfer  rate,  the  designer  should 
machine  as  many  fins  as  the  condenser  material  and  the  manu¬ 
facturing  process  will  allow.  A  saw  tooth  profile  is  pre¬ 
ferable  to  spacing  between  fins. 


I 


l 

I 

l 


DO  Form  1473 
S/J  01M  -014-6601 


2 


lieu***  CWM*l*l**Tt#**  T**'*  MMf*1* 


Approved  for  public  release;  distribution  unlimited 


Optimization  of  an  Internally 
Finned  Rotating  Heat  Pipe 


by 

William  A.  Davis,  Jr. 
Lieutenant  Commander,  United  States  Navy 
B.S.,  United  States  Naval  Academy,  1968 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  MECHANICAL  ENGINEERING 


from  the 

NAVAL  POSTGRADUATE  SCHOOL 
September  1980 


ABSTRACT 


A  finite  element  formulation  was  used  to  solve  the  steady- 
state  two-dimensional  conduction  heat  transfer  equation  in 
the  condenser  wall  section  of  an  internally  finned  rotating 
heat  pipe.  A  FORTRAN  program  using  this  method  was  coupled 
with  the  COPES/CONMIN  program  for  sensitivity  analysis  of 
design  variables  and  for  automated  design  of  the  internal 
heat  pipe  geometry. 

With  water  as  the  working  fluid,  numerical  results  ob¬ 
tained  for  copper  and  stainless  steel  heat  pipe  condenser 
sections  indicated  that  for  the  maximum  heat  transfer  rate, 
the  designer  should  machine  as  many  fins  as  the  condenser 

manufacturing  process  will  allow.  A  saw 
preferable  to  spacing  between  fins. 
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I.  INTRODUCTION 


A.  THE  ROTATING  HEAT  PIPE 

The  rotating  heat  pipe  is  a  closed  container  designed 
to  transfer  a  large  amount  of  heat  in  rotating  machinery. 

Its  three  main  component  parts  are:  a  cylindrical  evapora¬ 
tor,  a  truncated  cone  condenser,  and  a  working  fluid  as 
shown  in  Figure  1. 

At  rotation  above  the  critical  speed  of  a  rotating  heat 
pipe,  the  working  fluid  forms  an  annulus  in  the  evaporator, 
and  will  be  vaporized  by  heat  addition  to  it.  The  vapor 
flows  toward  the  condenser  as  a  result  of  a  pressure  differ¬ 
ence,  transporting  the  latent  heat  of  vaporization  with  it. 
External  cooling  of  the  condenser  causes  the  vapor  to  con¬ 
dense  on  the  inner  wall  and  release  its  latent  heat  of  eva¬ 
poration.  The  centrifugal  force  due  to  the  rotation  has  a 
component  acting  along  the  condenser  wall  that  will  act  to 
drive  the  condensate  back  to  the  evaporator  where  the  cycle 
is  repeated. 

In  a  conventional  heat  pipe,  the  force  driving  the  con- 
censate  back  to  the  evaporator  is  due  to  capillary  action, 
which  poses  a  limit  to  its  operation.  The  rotating  heat 
pipe  is  not  limited  by  capillary  action  and,  unlike  the 
thermosyphon  which  depends  on  gravity  to  cause  condensate 
return,  can  be  used  in  any  orientation  Cl]. 
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B.  OPERATING  LIMITS  OF  A  ROTATING  HEAT  PIPE 


The  first  theoretical  investigation  of  the  rotating 
heat  pipe  conducted  at  the  Naval  Postgraduate  School  was  per¬ 
formed  by  Ballback  [2]  in  1969-  He  studied  the  limitations 
of  performance  imposed  on  the  rotating  heat  pipe  due  to 
various  fluid  dynamic  mechanisms.  Using  existing  theory  and 
experimental  correlations,  he  was  able  to  estimate  the  sonic 
limit,  boiling  limit,  entrainment  limit,  and  the  condensing 
limit  of  performance. 

1 .  The  Sonic  Limit 

When  increasing  the  heat  flux  in  a  rotating  heat 
pipe,  it  is  possible  to  reach  a  limiting  flow  rate  of  the 
vapor  brought  on  by  a  choked  flow  condition  in  the  pipe. 

This  condition  imposes  a  limiting  value  on  the  amount  of 
energy  the  vapor  can  transport,  thus  reducing  the  effective¬ 
ness  of  the  heat  pipe.  The  limiting  heat  transfer  rate 
becomes 


<*t  =  pv  Uv  A  hfg 


(1) 


and  the  vapor  velocity  is  considered  to  be  sonic. 


uv  =  C  =  /  ggfcRT  (2) 

where 

Uv  =  velocity  of  the  vapor  in  ft/sec,  and 

2 

A  =  cross  sectional  area  for  the  vapor  flow  in  ft 
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c  =  sonic  velocity  in  ft/sec 

2 

gg  =  gravitational  constant,  32.1739  f t-lbm/lbf-sec 
k  =  ratio  of  specific  heats 
R  =  gas  constant  in  ft-lbf/lbm  ?.R,  and 
T  =  absolute  temperature  in  ?R. 

2 .  Boiling  Limit 

Kutateladze  [3]  postulated  that  the  transition  from 
nucleate  to  film  boiling  is  totally  a  hydrodynamic  process. 
He  determined  a  theoretical  formula  for  predicting  the  burn¬ 
out  flux 

Qt  =  K  ^  Ab  hfgta  g(Pf  -  PV>)1A  (3) 

where 

K  =  constant  value 

3 

Pv  =  density  of  the  vapor  in  lbm/ft 

2 

=  heat  transfer  area  in  the  boiler  in  ft 
h^g  =  latent  heat  of  vaporization  in  Btu/lbm 
o  =  surface  tension  in  lbf/ft 

2 

g  =  acceleration  of  gravity  in  ft/hr 

.  3 

p^.  =  density  of  fluid  m  lbm/ft 

.  3 

Pv  =  density  of  vapor  in  lbm/ft 

The  experimental  data  obtained  by  Kutateladze  suggested  a 

value  for  K  in  the  range  of  0-.  13  to  0.19. 
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3 .  Entrainment  Limit 


The  flooding  constraint  in  a  wickless  heat  pipe  was 
examined  by  Sakhuja  [4]  who  developed  the  correlation 


Q  .  Ax  C2  hfg/  g (4) 

17  (1  +  (P  V/Pf)l/4}2 

where 

Q  =  heat  transfer  rate  in  Btu/hr 

2 

=  flow  area  in  ft 

C  =  dimensionless  constant,  0.725  for  tube  with  sharp 
edged  flange 

h^  =  latent  heat  of  vaporization  in  Btu/lbm 

2 

g  =  acceleration  due  to  gravity  in  ft/hr 

D  =  inside  diameter  of  heat  pipe  in  ft 

3 

Pf  =  density  of  the  fluid  in  lbm/ft 

3 

Pv  =  density  of  the  vapor  in  lbm/ft 


4 .  Condensing  Limit 

Ballback  [2]  determined  the  condensation  solution 
for  a  rotating  heat  pipe  by  modeling  the  condenser  section 
of  a  rotating  heat  pipe  as  a  rotating  truncated  cone.  He 
developed  the  following  expression  for  the  condensation 
limit : 


Q 


t 


,2 

*{3 


■f  pf 


2  3  — 

oj  h  -  {T  -T  }  4 

fg  s  w 


sin  $ 


}  {CR 


8  8  3 

+  Lsin  0D3-Ro3}4 

(5) 
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where 


Q  =  total  heat  transfer  rate  in  Btu/hr 

kp  =  thermal  conductivity  of  the  condensate  film 

x  in  Btu/hr-f t-° F 

3 

p^.  =  density  of  fluid  m  lbm/ft 

u>  =  angular  velocity  in  1/hr 

h^  =  latent  heat  of  vaporization  in  Btu/lbm 

Ts  =  saturation  temperature  in  °F 

Tw  =  inside  wall  temperature  in  ?F 

y^  =  viscosity  of  fluid  in  Ibm/ft-hr 

^  =  half  cone  angle  in  degrees 

R  =  minimum  wall  radius  in  ft 

o 

L  =  length  along  the  wall  of  the  condenser  in  ft 

yf  =  viscosity  of  the  fluid  in  lbm/ft-sec 

The  condensing  limit  equation  is  a  function  of  the  geometry 
and  speed  of  the  rotating  heat  pipe,  and  the  physical  pro¬ 
perties  of  the  working  fluid. 

Tantrakul  [5]  calculated  these  limitations  for  a  heat 
pipe  with  specific  physical  characteristics  as  shown  in 
Table  I,  with  the  results  shown  in  Figure  2. 

TABLE  I 

Specification  of  a  Typical  Rotating  Heat  Pipe 


Length 

14.000 

inches 

Minimum  diameter 

2.000 

inches 

Wall  thickness 

0.125 

inches 

Internal  half  angle 

1.000 

degree 

Rotating  speed 

2700 

RPM 
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Obviously  from  the  results  in  Figure  2,  the  condensing  limit 

» 

is  the  predominant  limitation  for  the  amount  of  heat  that  can 
be  transferred  from  the  heat  pipe.  However,  the  other 
limitations  may  become  important  as  the  heat  pipe  geometry  and 
operating  conditions  are  varied. 

In  order  to  augment  the  heat  transfer  capacity  of  the 
heat  pipe,  recent  efforts  have  been  aimed  at  raising  the  con¬ 
densing  limit  line  which  may  be  accomplished  by: 

a.  a  high  value  of  cone  angle,  to  increase  the  centrifugal 
driving  force, 

b.  some  type  of  promoter  of  dropwise  condensation  to  in¬ 
crease  the  value  of  the  inside  heat  transfer  coeffi¬ 
cient,  h,  or 

c.  use  of  an  internally  finned  condenser  to  increase 

the  inner  wall  surface  area  and  the  value  of  h,  since 
the  presence  of  a  fin  will  decrease  the  effective 
condensate  film  thickness. 

A  high  value  of  cone  angle  means  a  departure  from  the 
cylinder  or  shaft  shape.  Since  the  principal  known  appli¬ 
cation  for  the  rotating  heat  pipe  is  in  the  cooling  cf  rota¬ 
ting  machinery  this  approach  was  not  pursued.  Although 
effective  promoters  of  dropwise  condensation  exist,  none  as 
yet  can  be  considered  to  be  permanent,  and  this  approach  was 
likewise  ruled  out.  The  remaining  alternative,  using  inter¬ 
nal  fins  in  the  condenser  section  to  raise  the  condensing 
limit,  was  seen  as  the  best  choice. 
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C.  ANALYSIS  OF  THE  INTERNALLY  FINNED  ROTATING  HEAT  PIPE 


Pursuing  the  addition  of  internal  fins  as  a  way  to  raise 
the  condensing  limit,  Schafer  [6]  developed  an  analytical 
model  for  a  heat  pipe  with  a  triangular  fin  profile  as  shown 
in  Figure  3.  He  assumed  one-dimensional  heat  conduction 
through  the  wall  and  fin.  Corley  [7]  for  this  same  case 
developed  a  two-dimensional  heat  conduction  model  using  a 
Finite  Element  Method,  and  also  assumed  a  parabolic  tempera¬ 
ture  distribution  along  the  fin  surface.  His  results  indi¬ 
cated  a  significant  improvement  in  heat  transfer  performance 
of  about  75%  above  that  predicted  by  the  one-dimensional  model 
of  Schafer  [6].  However,  Corley  [7]  cautiously  noted  a 
probable  error  of  50%  existed  at  the  fin  apex,  and  conse¬ 
quently  mentioned  that  there  may  be  a  total  heat  transfer 
error  of  as  high  as  15%.  Tantrakul  [5]  modified  Corley's 
computer  program  by  increasing  the  number  of  finite  elements 
in  order  to  minimize  the  heat  transfer  error  at  the  apex  of 
the  fin.  His  results  with  this  modification  converged  with 
the  results  of  Corley.  Purnomo  [1]  developed  a  two-dimen¬ 
sional  Finite  Element  Method  solution  using  a  linear  trian¬ 
gular  finite  element  model  as  shown  in  Figure  4.  Purnomo ' s 
[1]  Finite  Element  Method  program  also  worked  and  converged. 
Purnomo ' s  [1]  code,  when  made  to  approach  the  geometry  of  a 
smooth  tube,  did  not  agree  with  the  analytical  and  experi¬ 
mental  data  obtained  by  Schafer  [5]  for  a  smooth  tube.  This 
cast  doubts  about  the  validity  of  Purnomo ' s  code. 


The  parametric  studies  conducted  using  Purnomo ' s  code 
gave  no  clear  indication  of  the  best  condenser  geometry  to 
maximize  heat  transfer.  Also,  his  code  was  tedicus  to  use 
and  required  numerous  runs  to  obtain  data  since  it  was 
written  to  perform  only  one  analysis  at  a  time. 

Clearly  a  computer  program  that  could  make  numerous  runs 
with  minimal  data  input  and  could  also  automatically  find 
improved  designs  would  be  valuable. 


D.  THESIS  OBJECTIVES 

The  objectives  of  this  thesis  were  therefore: 

1.  To  modify  Purnomo ' s  Cl]  computer  program  so  that  it 
is  compatible  with  the  COPES/COMNIN  program  [S]  and 
can  be  used  for  analysis  and  automated  design  of 
rotating  heat  pipes  (internally  finned  or  smooth). 

2.  To  compare  results  using  Purnomo ' s  code  with  analyti¬ 
cal  results  for  a  smooth  tube  obtained  by  Schafer  [5] 
to  determine  if  and  where  an  error  exists. 

3.  To  use  the  sensitivity  analysis  capability  of  COPES/ 
CONMIN  in  conjunction  with  the  modified  program  to 
study  heat  transfer  in  an  internally  finned  rotating 
heat  pipe. 

4.  To  use  the  resulting  program  to  obtain  an  optimum 
design  for  an  internally  finned  rotating  heat  pipe 
to  obtain  experimental  data  to  compare  with  the 
analytical  results. 

5.  To  use  the  resulting  program  to  obtain  numerical 
results  in  place  of  data  obtained  from  expensive 
experimental  operations. 
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II.  NUMERICAL  OPTIMIZATION 


A.  BACKGROUND 

Most  design  processes  require  the  minimization  or  maxi¬ 
mization  of  some  parameter  which  may  be  called  the  design 
objective.  For  the  design  to  be  acceptable,  it  must  satisfy 
a  variety  of  physical,  aesthetic,  economic  and,  on  occasion, 
political  limitations  which  are  referred  to  here  as  design 
constraints.  While  part  of  the  design  problem  may  not  be 
easily  quantified,  most  of  the  design  criteria  can  be 
described  in  numerical  terms. 

To  the  extent  that  the  problem  can  be  stated  in  numerical 
terms,  a  computer  program  can  be  written  to  perform  the 
necessary  calculations.  For  this  reason,  computer  analysis 
is  commonplace  in  most  engineering  organizations.  For 
example,  in  structural  design  the  configuration,  materials, 
and  loads  may  be  defined  and  a  finite  element  analysis  com¬ 
puter  code  is  used  to  calculate  stresses,  deflections,  and 
other  response  quantities  of  interest.  If  any  of  these 
parameters  are  not  within  prescribed  bounds,  the  engineer 
may  change  the  structural  member  sizes  and  rerun  the  program. 
The  computer  code  therefore  provides  only  the  analysis  of  a 
proposed  design,  with  the  engineer  making  the  actual  design 
decisions.  This  approach  to  design,  which  may  be  called 
computer-aided  design,  is  commonly  used  today. 


Another  common  use  of  analysis  codes  is  in  tradeoff 
studies.  For  example,  an  aircraft  trajectory  analysis  code 
may  be  run  repetitively  for  several  payloads,  calculating 
the  aircraft  range,  to  dfetermine  the  range-payload  sensitivity 

A  logical  extension  to  computer-aided  design  is  fully 
automated  design,  where  the  computer  makes  the  actual  design 
decisions,  or  performs  trade-off  studies  with  a  minimum  of 
man-machine  interaction  [9]. 

B.  CONSTRAINED  FUNCTION  MINIMIZATION  (CONMIN) 

Vanderplaats  [10]  developed  an  optimization  program  CCNMIN 
capable  of  optimizing  a  very  wide  class  of  engineering  pro¬ 
blems.  CONMIN  is  a  FORTRAN  program,  in  subroutine  form,  that 
optimizes  a  multi-variable  function  subject  to  a  set  of  in¬ 
equality  constraints  based  on  Zoutendijk’s  [11]  method  of 
feasible  directions  [12]. 

Three  basic  definitions  are  required  to  discuss  the 
use  of  CONMIN: 

Design  Variables  -  Those  parameters  which  the  optimization 
program  is  permitted  to  change  in  order  to  improve  the 
design.  Design  variables  appear  only  on  the  right  hand 
side  on  an  equation  and  are  continuous. 

Design  Constraints  -  Any  parameter  which  must  not  exceed 
specified  bounds  for  the  design  to  be  acceptable. 

Design  constraints  may  be  linear  or  nonlinear,  implicit 
or  explicit,  but  they  must  be  continuous  functions  of 
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the  design  variables.  Design  constraints  appear  only 
on  the  left  side  of  the  equations. 

Objective  Function  -  The  parameter  which  is  going  to  be 
minimized  or  maximized  during  the  optimization  process. 
The  objective  function  may  be  linear  or  nonlinear, 
implicit  or  explicit,  and  must  be  a  continuous  function 
of  the  design  variables.  The  objective  function  usually 
appears  on  the  left  side  of  an  equation,  but  it  may 
appear  on  the  right  side  if  it  is  also  a  design  variable. 
Design  constraints  and  objective  functions  are  usually  inter¬ 
changeable  . 

C.  CONTROL  PROGRAM  FOR  ENGINEERING  SYNTHESIS  (COPES) 

Recall  that  the  optimization  program,  CONMIN,  was  written 

in  subroutine  form.  Vanderplaats  [8]  has  developed  a  main 
program  to  simplify  the  use  of  CONMIN  and  to  further  aid  in 
the  optimization  process.  The  user  must  supply  an  analysis 
subroutine  named  ANALIZ.  What  follows  are  programming  guide¬ 
lines  to  ensure  compatability  with  COPES. 

D.  PROGRAMMING  GUIDELINES 

In  developing  any  computer  code  for  engineering  analysis, 
it  is  prudent  to  write  the  code  in  such  a  way  that  it  is 
easily  coupled  to  a  general  synthesis  program  such  as  COPES. 
Therefore,  a  general  programming  practice  is  outlined  here 
which  in  no  way  inhibits  the  use  of  the  computer  program  in 
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its  traditional  role  as  an  analytic  tool,  but  allows  for 
simple  adaption  to  COPES.  This  approach  is  considered  good 
programming  practice  and  provides  considerable  flexibility 
of  design  options.  Only  five  basic  rules  must  be  followed:. 

I.  Write  the  code  in  subroutine  form  with  the  primary 
routine  called  as  SUBROUTINE  ANALIZ ( ICALC) .  The 
name  ANALIZ  is  compatable  with  the  COPES  program  and 
ICALC  is  a  calculation  control.  Note  that  subroutine 
ANALIZ  may  call  numerous  other  subroutines  as  required 
to  perform  the  necessary  calculations. 

II.  Segment  the  program  into  INPUT,  EXECUTION,  and  OUTPUT. 
The  calculation  control,  ICALC,  will  determine  the 
portion  of  the  analysis  code  to  be  executed.  ICALC=1; 
the  program  reads  all  data  required  to  perform  the 
analysis.  Also,  any  initialization  of  constants 
which  will  be  used  repetively  during  execution  is 
done  here.  This  initial  input  information  is  printed 
here  for  later  reference  and  for  program  debugging. 
ICALC=2;  the  program  performs  the  execution  phase  of 
the  analysis  task.  No  data  reading  or  printing  is 
done  here,  except  on  user-defined  scratch  disc.  Data 
may  be  printed  here  during  program  debugging,  in 
which  case  it  should  be  controlled  by  a  print  con¬ 
trol  parameter  which  is  read  during  input.  In  this 
way,  this  print  may  be  turned  off  after  the  program 
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is  debugged,  but  may  be  used  again  during  future 
program  expansion  debugging.  The  reason  that  print¬ 
ing  is  not  allowed  during  execution  is  that  when 
optimization  is  being  done,  the  code  will  be  called 
many  times  with  ICALC=2,  resulting  in  voluminous 
print.  ICALC=3;  the  results  of  the  anlaysis  are 
printed.  Also  the  essential  input  parameters  which 
may  have  been  changed  during  optimization  should  be 
printed  here  for  easy  reference.  In  summary,  when: 
ICALC  =  1  Read  input  data. 

ICALC  =  2  Execute  the  analysis. 

ICALC  =  3  Print  the  results. 

Store  all  parameters  which  may  be  design  variables, 
objective  functions  or  constraints  in  a  single 
labeled  common  block  called  GLOBCM.  The  order  in 
which  they  are  stored  is  arbitrary.  A  listing  of  the 
COPES  program  should  be  checked  to  see  how  many  para¬ 
meters  may  be  stored  in  GLOBCM  (the  dimension  of 
ARRAY) .  Initial  distribution  of  COPES  allows  for  1500 
parameters . 

During  execution  or  output,  no  parameters  which  are 
read  during  input  should  be  updated.  For  example, 
if  variable  X  is  initialized  during  input,  the 
execution  segment  must  not  update  X  such  as  X=X  + 

3.2.  Instead  a  new  variable,  Y=X  +  3.2  must  be 
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V.  Write  all  programs  in  standard  language,  avoiding 
machine  dependent  capabilities  such  as  seven  letter 
FORTRAN  names  (CDC).  While  this  guideline  is  not 
essential  to  the  use  of  the  analysis  code  within  the 
COPES  program,  it  makes  the  analysis  code  much  more 
transportable  between  different  computer  systems,  a 
capability  which  easily  justifies  a  slight  reduction 
in  efficiency  on  a  given  machine. 

Adherence  to  these  guidelines  not  only  leads  to  a  more 
readable  and  machine  independent  computer  code,  but  allows 
this  code  to  be  coupled  to  the  COPES  program  without  modifi¬ 
cation. 

Having  written  the  analysis  code,  it  may  be  executed 
either  with  a  simple  main  program  or  within  the  COPES  pro¬ 
gram  to  perform  the  analysis.  To  insure  that  guideline  IV 
is  followed,  the  following  main  test  program  is  recommended. 
Note  that  this  program  calls  ANALIZ  twice  with  ICALC=2  and 
ICALC=3,  to  show  that  the  same  result  is  obtained  repeti¬ 
tively  . 

C  MAIN  PROGRAM  TO  CHECK  SUBROUTINE  ANALIZ. 

C  READ,  EXECUTE ,  AND  PRINT  ^ 

DO  10  ICALC=1 , 3 
10  CALL  ANALIZ  (ICALC) 

C  EXECUTE  AND  PRINT  AGAIN  TO  BE  SURE  THE  RESULTS 
C  DO  NOT  CHANGE 
DO  20  ICALC  =  2 , 3 
20  CALL  ANALIZ  (ICALC) 

STOP 

END 


III.  FINITE  ELEMENT  SOLUTION 


A.  REVIEW  OF  THE  PREVIOUS  ANALYSIS 

Schafer  [6]  studied  the  one-dimensional  model  heat  trans¬ 
fer  solution  and  Corley  [7]  studied  the  two-dimensional 
model  for  an  internally  finned  rotating  heat  pipe.  Both 
used  the  same  assumptions  and  boundary  conditions  based  upon 
the  analysis  of  3allback  [2],  which  are  similar  to  these  used 
in  the  Nusselt  analysis  of  film  condensation  on  a  flat  wall. 

The  more  important  of  those  assumptions  are: 

1.  steady  state  operation, 

2.  film  condensation,  as  opposed  to  dropwise  condensation, 

3 .  laminar  flow  of  the  condensate  film  along  both  the 
fin  and  the  trough, 

4.  static  balance  of  forces  within  the  condensate, 

5.  one -dimensional  conduction  heat  transfer  through 
the  film  thickness  (no  convective  heat  transfer  in 
the  condensate  film), 

6.  no  liquid  -  vapor  interfacial  shear  forces, 

7.  no  condensate  subcooling, 

8.  zero  heat  flux  boundary  conditions  on  both  sides  of 
the  wall  section  (symmetry  conditions),  as  shown 

in  Figure  5, 

9.  saturation  temperature  at  the  fin  apex, 
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10.  zero  film  thickness  at  the  fin  apex,  and 

11.  negligible  curvature  of  the  condenser  wall. 

Purnomo  [1]  developed  a  two-dimensional  Finite  Element 

Method  solution  using  a  linear  triangular  finite  element 
model  as  shown  in  Figure  4.  Purnomo  modif ied  Corley's  assump¬ 
tion  that  the  fin  apex  was  at  the  saturation  temperature  and 
allowed  the  value  of  the  temperature  at  the  apex  to  float. 

He  assumed  a  parabolic  temperature  distribution  along  the 
fin  surface. 

Purnomo ' s  statement  of  the  problem  for  the  formulation 
of  the  Finite  Element  Method  as  shown  in  Figure  5  is 


2  2 

a  t  +  <rr  .  Q 


ax"  ay 


(6) 


with  the  boundary  conditions: 


3T  _ 


a)  along  boundary  1,  -k  ^  =  h^  (T  -  T 


a  t 


b)  along  boundary  3,  -k  j-pj-  =  h2  "  Tco^ 


c)  along  boundaries  2  and  4,  5— 

d  n 


a  t  . 


A  detailed  description  of  the  numerical  formulation  is  pre¬ 
sented  in  his  thesis. 
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Purnomo' s  computer  program  consisted  of  a  main  program 
and  three  subroutines; 

a)  the  routine  "COORD"  used  to  define  positions  of 
system  coordinate  points, 

b)  the  routine  "FORMAT"  used  to  formulate  the 
Finite  Element  Method  equations,  and 

c)  the  routine  "BANDEC"  as  an  equation  solver  for  a 
symmetric  matrix  which  has  been  transformed  into 
banded  form. 


B.  THE  COMPUTER  PROGRAM  DEVELOPMENT 

Purnomo ' s  [1]  two-dimensional  finite  element  program 
is  the  basis  for  the  present  analysis  code.  The  first 
task  undertaken  in  the  development  of  this  thesis  was  to 
check  Purnomo 's  code  for  validity.  Theoretical  results 
using  Purnomo 's  code  for  the  case  of  zero  fins  were  compared 
to  theoretical  results  and  experimental  data  for  a  smooth 
tube  obtained  by  Schafer  [6].  Purnomo ' s  results  exceeded 
Schafer's  results  by  a  factor  of  approximately  two.  In 
studying  this  discrepancy,  an  error  was  discovered  in 
Purnomo 's  code.  In  encoding  formula  (11-11)  [1]  for  mass 
flow  rate  in  FORTRAN,  the'  sin  <j>  ’  term  was  dropped.  The 
correct  form  of  the  equation  is  shown  below. 


.  p 2 ( Rn  +  xsin  <J>)<$*(x)2  sin  <j>  _ 

M  (x)  =  - - - [5  * ( x ) £  +5  * ( X )  ta n  a] 

x  3  uf 
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The  effect  of  this  error  was  that  the  film  thickness  along 
the  length  of  the  condenser  remained  very  small.  Purnomo ' s 
original  code  consequently  indicated  a  much  higher  value  of 
heat  transfer  rate  than  the  correct  analysis  did.  In  a 
correct  analysis  code,  the  condensate  film  thickness  grows 
continuously  until  the  condensate  reaches  the  evaporator.  A 
plot  of  the  condensate  film  thickness  calculations  from 
Purnomo ' s  original  code  for  the  case  of  a  smooth  tube,  his 
corrected  code,  and  Schafer's  program  for  a  smooth  tube  are 
shown  in  Figure  6.  Purnomo ' s  corrected  code  for  the  case  of 
zero  fins  agreed  to  within  8  percent  of  the  results  for 
heat  transfer  rate  using  Schafer's  [6]  analysis  for  a  smooth 
tube . 

Since  no  experimental  data  existed  for  further  comparison 
of  the  finned  model,  the  next  task  undertaken  was  to  adapt 
the  analysis  code  to  permit  automated  design  and  sensitivity 
analysis  using  COPES/CONMIN .  Many  modifications  were  made, 
some  of  which  are  mentioned  here.  The  program  was  rewritten 
in  subroutine  form  and  segmented  into  input,  execution,  and 
output  sections  to  make  it  compatible  with  COPES.  Since 
COPES  was  written  to  use  single  precision  mathematics  and 
the  analysis  code  uses  double  precision  to  allow  for  possible 
ill  conditioning,  the  subroutine  ANALIZ  also  makes  the  trans¬ 
formation  from  single  to  double  precision.  The  initial  value 
of  the  film  thickness  is  calculated  within  the  code  based  on 
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a  formulation  by  Sparrow  and  Gregg  [13].  Previously  a 
constant  value  based  on  the  same  analysis  was  used. 

The  modified  code  can  be  used  alone  for  analysis  of  a 
given  geometry  or  it  can  be  used  in  conjunction  with  COPES/ 
CONMIN  for  a  single  analysis,  sensitivity  analysis,  or 
automated  design  (optimization) . 

A  listing  of  the  revised  computer  program  is  included 
as  Appendix  B. 


C.  SENSITIVITY  ANALYSIS 

The  sensitivity  analysis  feature  of  COPES/CONMIN  was  used 
to  obtain  data  for  the  heat  transfer  rate  Q  as  a  function 
of  various  design  variables.  All  design  variables  were  held 
constant  except  the  one  of  interest  which  was  varied  as 
specified  in  block  Q  of  the  COPES  input  data. 

For  example,  to  obtain  the  heat  transfer  rate  Q  for  26 
different  values  of  rotational  speed,  the  following  data 
input  is  required: 


$  31ock  P 
1 

15 

$  Block  Q 
16,  26 

3000.,  100.,  200.,  300.,  400.,  500.,  600.,  700., 

300.,  900.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000., 

7000.,  3000.,  9000.,  10000.,  11000.,  12000.,  13000., 

14000.,  15000.,  3000. 
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For  each  plot  only  the  variable  on  the  ordinate  is 
changed.  All  other  design  variables  remain  constant.  The 
basic  design  used  for  sensitivity  plots  has  the  design 
variable  values  shown  below: 


TABLE  II 

BASIC  DESIGN  FOR  SENSITIVITY  ANALYSIS 


condenser  length 

CLI 

= 

9.0  inches 

cone  half  angle 

CANGL 

= 

1.0  DEGREES 

condenser  radius 

RBASEI 

= 

0.775  inches 

wall  thickness 

THICKI 

= 

0.03125  inches 

fin  height 

B 

= 

0.025  inches 

speed  of  rotation 

RPM 

= 

3000  RPM 

saturation  temperature 

TSAT 

r 

100°  F 

ambient  temperature 

TINF 

= 

60°  F 

outside  heat  transfer 
coefficient 

HINF 

= 

5000  Btu/HR  FT 

fin  half  angle 

FANGL 

= 

10 

ratio  of  trough  width 
to  fin  base  width 

ZOA 

r 

12.8 

number  of  fins 

ZFIN 

= 

40.0 

Note:  In  the  code,  ZFIN  appears  on  the  left  hand  side  of 

an  equation  and  is  therefore  by  strict  definition  not  a 
"design  variable".  ZFIN  is  calculated  from  the  values  of 
CBASE ,  EZERO,  and  EPSO. 

D.  DESIGN  OPTIMIZATION 

Counting  parameters  such  as  external  heat  transfer 
coefficient  HINF,  there  are  thirteen  possible  design 
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variables.  Nine  of  these  are  geometric  or  functional  para¬ 
meters  such  as  wall  thickness,  fin  height,  and  speed  of 
rotation.  The  design  variables,  possible  constraint  func¬ 
tions,  and  the  objective  function  appear  in  the  Global 
common  block  GL0B1  in  the  code  and  are  listed  below  by  fortran 
name  for  clarity. 

DESIGN  VARIABLES 
CLI 
CANGL 
RBASEI 
R2I 

THICKI 

BFIN 

TZ 

TSS 

TINF 

HINF 

FANGL 

ZOA 

RPM 

CONSTRAINT  FUNCTIONS 
ZFIN 
BOA 
DIFF 

OBJECTIVE  FUNCTION 
QTOT 


There  are  a  wide  variety  of  design  problems  that  can  be 
pursued  with  the  code.  For  example  one  might  wish  to 
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determine  the  smallest  length  condenser  section  and  the  best 
internal  geometry  for  a  specified  heat  transfer  rate-- 
perhaps  to  cool  an  electric  motor. 

To  maximize  heat  transfer  rate  through  the  condenser 
wall  the  designer  adds  a  number  of  fins  to  increase  the 
inside  surface  area.  As  more  and  more  fins  are  added 
however,  the  cross-sectional  area  for  conduction  through 
each  fin  is  decreased.  Also,  the  film  thickness  of  the 
condensate  in  the  though  increases  and  in  fact  could  com¬ 
pletely  cover  the  fins  and  substantially  reduce  heat  trans¬ 
fer  through  the  fin.  So  there  should  exist  some  optimum 
combination  of  number  of  fins,  fin  height,  fin  half  angle, 
and  ratio  of  trough  width  to  fin  width  that  will  permit 
maximum  heat  transfer  rate. 

The  design  study  undertaken  was  to  determine  the  fin 
height,  fin  half  angle,  and  fin  spacing  which  would  yield 
the  maximum  heat  transfer  possible.  The  design  variables 
then  were  BFIN,  FANGL,  and  ZOA.  Other  potential  design 
variables  were  held  constant.  The  objective  function  to  be 
maximized  was  QTOT ,  heat  transfer  rate  out  of  the  condenser. 

For  comparison,  the  theoretical  upper  limit  on  heat 
transfer  was  calculated  based  on  an  external  surface  tempera¬ 
ture  equal  to  the  working  fluid  saturation  temperature.  This 
assumes  that  there  is  no  thermal  resistance  across  the 
condensate  and  the  condenser  wall.  When  this  upper  limit 
was  used,  the  maximum  heat  transfer  rate  was  predicted  to  be 
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63,322  BTU/HR  using  the  following  formula: 


Q  a  h  2 nr  L  (Twall  -  T«)  (9) 

max 

where 

h  =  outside  convective  heat  transfer  coefficient 
(5000  BTU/HR-FT2-°F) 

r  =  average  outside  radius  of  condenser  wall  (0.0056  ft) 
L  =  condenser  length  (0.75  feet) 

Twall=temperature  of  the  outside  wall  (100°F) 

T«  =  ambient  temperature  (60°F) 

Certain  constraints  were  placed  on  the  design  based  on 
engineering  judgement.  For  example,  the  number  of  fins  was 
not  allowed  to  exceed  400  and  the  minimum  fin  half  angle 
allowed  was  10  degrees.  These  values  were  based  on  structural 
and  manufacturing  considerations. 


IV.  RESULTS 


A.  INTRODUCTION 

The  computer  code  was  used  in  conjunction  with  the  COPES/ 
CONMIN  program  for  sensitivity  analysis  and  design  optimiza¬ 
tion  for  maximum  heat  transfer  rate.  Numerical  results  are 
discussed  below. 

B.  SENSITIVITY  ANALYSIS 

As  shown  in  Figures  7  and  8,  the  heat  transfer  rate  in¬ 
creases  for  an  independent  increase  of  the  design  variable 
of  interest,  but  levels  off  at  the  theoretical  maximum  heat 
transfer  rate  for  the  given  overall  geometry.  The  heat  trans¬ 
fer  rate  then  is  limited  by  the  external  resistance  which  has 
become  the  controlling  factor. 

In  Figure  9,  the  heat  transfer  rate  increases  in  a  simi¬ 
lar  way  for  an  increase  in  external  heat  transfer  coefficient. 
Again,  the  rate  of  increase  appears  to  lessen  for  external 
heat  transfer  coefficients  above  10,000  BTU/HR- FT^ • ° F  due  to 
other  limiting  resistances.  In  Figure  10  the  heat  transfer 
rate  is  observed  to  increase  linearly  with  condenser  length. 

This  is  expected  since  0  is  a  function  of  the  area  and  for  small 
values  of  t>  the  area  varies  directly  with  length.  In  Figure  11, 
Q  is  seen  to  rise  in  a  non-linear  manner  for  an  increase  in 


37 


cone  half  angle,  <p .  This  is  because  the  internal  heat  trans¬ 
fer  coefficient,  h^  is  a  non-linear  function  of  cone  half 
angle . 

Purnomo  [1]  concluded  that  the  heat  transfer  rate  continu¬ 
ously  increased  as  the  fin  half  angle  decreased  and  that  this 
was  largely  a  result  of  the  fact  that  the  number  of  fins  in¬ 
creased  at  the  same  time.  He  also  stated  that  the  increase 
in  heat  transfer  was  only  slight  when  the  fin  half  angle  was 
less  than  11  degrees.  This  is  somewhat  misleading  in  that  in 
his  analysis  the  number  of  fins  was  being  changed  with  every 
change  in  fin  half  angle.  This  author  has  drawn  a  different 
conclusion.  Figure  12  shows  a  plot  of  heat  transfer  rate  vs. 
varying  fin  half  angle  all  for  a  condenser  with  40  fins.  The 
heat  transfer  rate  as  a  function  of  fin  half  angle  rises 
sharply  from  1-11  degrees  and  continues  to  rise,  but  less 
steeply,  as  the  fin  half  angle  increases  beyond  11  degrees. 

The  maximum  fin  half  angle  possible  with  40  fins  is  67.5 
degrees . 

C.  CONSTRAINED  OPTIMIZATION 

In  the  design  problem  undertaken  to  determine  the  optimum 
internal  geometry  for  maximum  heat  transfer,  numerous  runs 
were  made  for  condensers  made  of  copper,  stainless  steel,  and 
a  ceramic  material.  These  materials  have  thermal  conduc¬ 
tivity  values  of  231,9,  and  1.0  BTU/HR*FT’°F  respectively. 

It  was  expected  that  for  each  material  a  different  optimum 
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design  would  emerge.  The  results,  therefore,  were  unexpected. 
From  different  starting  points  (original  designs),  some  out¬ 
side  the  feasible  region,  and  for  external  heat  transfer 

2 

coefficients  from  1000-50,000  BTU/HR • FT" • 0 F ,  the  same  optimum 
design  for  maximum  heat  transfer  was  reached,  which  is  out¬ 
lined  in  Table  III  below.  Each  material,  of  course,  had  a 
different  heat  transfer  rate  even  though  the  geometry  for 
maximum  heat  transfer  was  the  same. 

TABLE  III 

CONSTRAINED  OPTIMIZATION  RESULTS  FOR  COPPER 
STAINLESS  STEEL,  AND  CERAMIC  MATERIAL 

Optimum  Design  for  All  Materials 
fin  height  0.023  inches 


fin  half  angle 

10.0  degrees 

number  of  fins 

400 . 0 

ratio  of  trough  width 

0.5 

to  fin  base  width 

MATERIAL 

^external 
( BTU/HR • FT2 • 0  F) 

HEAT  TRANSFER  RATE 
(BTU/HR) 

Copper 

1000 

13,822 

50,000 

296,560 

Stainless  Steel 

1000 

10,500 

50,000 

38,650 

Ceramic  Material 

1000 

3844 

50,000 

52  37 

D.  UNCONSTRAINED  OPTIMIZATION 


The  results  of  the  constrained  optimization  runs  were  so 
unexpected  that  it  was  decided  to  remove  the  constraint  on 
the  humber  of  fins  and  the  trough  width  to  fin  width  ratio 
and  repeat  the  optimization  runs.  If  another  identical  design 
was  again  reached  for  all  materials,  the  code  would  have  to  be 
considered  in  error.  The  results  are  presented  in  Table  IV. 

The  optimum  design  for  maximum  heat  transfer  was  different 
for  each  material  and  each  heat  transfer  coefficient.  For 
each  material,  the  optimum  fin  height  b  increases  as  the  ex¬ 
ternal  heat  transfer  coefficient  increases.  This  provides 
more  fin  surface  area  to  increase  heat  transfer  rate  as  the 
outside  heat  transfer  resistance  no  longer  controls.  Despite 
a  fifty-fold  change  in  heat  transfer  coefficient,  each 
material  maintained  essentially  a  constant  5 */b  ratio.  The 
higher  conductivity  materials  required  less  exposed  fin  sur¬ 
face  than  the  lower  conductivity  materials.  Copper,  for 
instance,  could  be  57  percent  covered  by  the  trough  condensate 
while  the  ceramic  material  was  only  14-19  percent  covered  for 
maximum  heat  transfer. 

In  each  material  the  number  of  fins  decreased  for  an 
increase  in  external  heat  transfer  coefficient.  This  provides 
more  space  in  the  troughs  to  carry  the  increased  condensate 
which  results.  Otherwise,  fin  performance  would  be  degraded. 
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UNCONSTRAINED  OPTIMIZATION  RESULTS 
FOR  COPPER,  STAINLESS  STEEL,  AND  CERAMIC  MATERIAL 
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For  an  increase  in  external  heat  transfer  coefficient  the 
highest  percentage  change  in  heat  transfer  rate  occurs  in  the 
material  with  the  highest  thermal  conductivity.  In  the 
ceramic  material,  for  instance,  a  fifty-fold  increase  in  ex¬ 
ternal  heat  transfer  coefficient  results  in  only  a  36  percent 
increase  in  heat  transfer  rate.  In  the  copper  condenser,  the 
heat  transfer  rate  is  increased  by  a  factor  of  21.  This  shows 
the  strong  dominance  of  the  wall  resistance  to  heat  transfer 
in  the  ceramic  material. 

For  all  three  materials,  the  trough  is  essentially  eli¬ 
minated  for  the  optimum  design.  Why  the  trough  is  not 
entirely  eliminated  in  the  stainless  steel  condenser  section 
is  not  clear. 

Because  the  manufacture  of  a  very  large  number  of  fins  in 
a  small  diameter  condenser  is  not  practical,  it  is  desirable 
to  use  the  constrained  optimization  results  as  a  design  guide¬ 
line  . 

A  comparison  of  the  unconstrained  and  constrained  design 
for  a  copper  condenser  section,  for  instance,  shows  that 
using  the  more  realistic  number  of  400  fins  and  its  proper 
fin  height  b  instead  of  2040  fins  results  in  a  heat  transfer 
rate  difference  of  only  0.6  percent.  This  is  shown  clearly 
in  Table  V. 
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TABLE  V 

Comparison  of  Unconstrained  and  Constrained  Optimum  Design  for 
a  Copper  Condenser  Section  with  External^Heat  Transfer  Coeffi¬ 
cient  of  50000  BTU/HR-FT  -°F. 


constrained 


b  (inches) 
0.023 


Number  %Difference 

of  Fins  Q  (BTU/HR)  in  Q 

4Q0  296,560  - 


unconstrained 


0.019 


2040  298,330  0.6% 


The  author's  conclusion  is  that  within  the  realm  of  what 
can  now  be  manufactured,  the  same  basic  design  is  best  for 
all  materials  regardless  of  the  external  heat  transfer 
coefficient.  That  is,  the  designer  should  machine  as  many 
fins  as  the  condenser  material  and  the  manufacturing  process 
will  allow.  The  code  should  be  used  to  determine  the  fin 
height  to  avoid  degrading  the  fin  efficiency  by  too  high  a 
level  of  condensate  in  the  trough.  As  seen  in  Table  V,  it 
may  be  possible  to  lower  the  constraint  on  the  number  of  fins 
more  than  once  and  compare  the  resulting  heat  transfer  rate 
to  achieve  an  effective  but  less  expensive  design  to  manu¬ 
facture. 

E.  A  CAUTIONARY  NOTE 

The  reader  should  be  cautioned  that  when  an  analysis  code 
based  on  assumptions  made  for  certain  conditions  is  linked  to 
an  optimizer,  the  code  may  change  the  geometry  or  other 
conditions  such  that  the  original  assumptions  are  no  longer 
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valid.  Judicious  use  of  side  constraints  and  engineering 
common  sense  can  prevent  disaster.  Numerical  solutions  must 
not  be  blindly  accepted. 

It  should  be  noted  that  all  calculations  in  this  thesis 
were  made  with  water  as  the  working  fluid.  Sensitivity  data 
were  run  with  an  original  rotational  speed  of  3C00  RPM  while 
design  optimization  was  done  for  a  rotational  speed  of  3600 
RPM.  Different  combinations  of  operating  parameters,  working 
fluids,  and  condenser  section  materials  may  not  behave  in  a 
predictable  manner.  While  certain  trends  are  shown  for  the 
particular  condenser  section  studied  here ,  it  would  be 
preferable  to  modify  the  code  to  analyze  the  particular 
problem  at  hand  than  to  extrapolate  these  findings  to  too 
broad  an  application.  ' 

While  a  design  improvement  is  almost  certain  using  an 
optimization  routine  like  COPES/CONMIN ,  there  is  no  guarantee 
that  a  global  optimum  has  been  found.  Therefore,  for  a  given 
design  problem,  the  designer  should  start  the  optimization 
process  from  several  different  original  designs  and  compare 
the  results. 
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V.  CONCLUSIONS 


1.  The  computer  code  is  valid  and  can  be  used  for  single 
analysis,  sensitivity  analysis,  and  automated  design  of  an 
internally  finned  rotating  heat  pipe.  It  can  also  be  used  to 
generate  numerical  data  to  find  a  correlation  useful  for 
design. 

2.  For  zero  fins,  the  code  converges  to  the  results  for 
a  smooth  tube  obtained  by  Schafer  [6]. 

3.  The  heat  transfer  rate  increases  for  an  independent 
increase  in  fin  half  angle,  rotational  speed,  and  number  of 
fins  but  levels  off  at  the  theoretical  maximum  heat  transfer 
rate  for  the  heat  pipe.  The  heat  transfer  rate  continues  to 
increase  for  an  increase  in  condenser  length  and  cone  half 
angle . 

4.  Within  the  realm  of  known  materials,  maximum  heat 
transfer  occurs  for  the  same  fin  geometry  regardless  of  the 
external  heat  transfer  coefficient.  That  is,  for  a  given 
condenser  radius,  one  should  machine  as  many  fins  as  possible 
with  the  condenser  material  used  to  maximize  heat  transfer. 
The  code  should  be  used  to  determine  the  correct  fin  height 
to  avoid  degradation  of  fin  performance  by  too  high  a  level 
of  condensate  in  the  trough. 
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5 .  When  the  constraint  on  minimum  fin  half  angle  is  re¬ 
moved  the  optimum  design  is  a  large  number  of  very  thin  fins. 
A  practical  design  then  might  be  a  series  of  long  rectangular 
fins . 
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VI.  RECOMMENDATIONS 


Based  upon  the  calculated  results  of  this  thesis,  the 
following  recommendations  are  made. 

1.  Build  an  internally  finned  heat  pipe  to  obtain 
experimental  data  for  comparison  to  the  analytical 
predictions. 

2.  Generate  a  new  model  to  analyze  the  case  where  cone 
half  angle  is  zero.  This  would  greatly  decrease  the 
cost  of  manufacture,  and  permit  the  use  of  internally 
finned  heat  exchanger  tubing  which  is  commercially 
available . 

3.  Analyze  different  shaped  fins,  including  rectangular. 

4.  Obtain  numerical  data  using  the  code  to  attempt  to 
find  a  correlation  for  the  heat  transfer  rate  as  a 
function  of  important  variables: 

Q  =  Q  (temperature,  condenser  geometry,  fluid 

properties,  condenser  material  properties, 
heat  transfer  coefficients,  etc.) 

5.  Modify  the  code  to  include  as  constraint  functions 
structural  failure  modes  such  as  rupture  due  to  the 
internal  spinning  mass,  buckling,  and  whirling. 
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Figure  1.  Schematic  Drawing  of  a  Rotating  Heat  Pip 


SONIC  LIMIT 


ENTRAINMENT  LIMIT 


SO  180  200 

fs.nperature  Ts  deg  F 


Figure  3.  Internally  Finned  Condenser  Geometry, 
Showing  Fins,  Troughs,  and  Lines  of 
Symmetry 
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Figure  4.  Condenser  Geometry  Considered  with  25 
Linear  Triangular  Finite  Elements 


T 

i!l  +  III 

3x2  3y2 


b.c . 

a)  -  k  (T  -  Tgat)  Along  Boundary  Cl] 

b)  -  k  *  h2  (T  -  T^)  Along  Boundary  [3] 

3  T 

c)  r —  =0  Along  Boundaries  [2]  and  [4] 

3  n 

Figure  5.  Differential  Equation  and  Boundary  Conditions 
Considered  in  the  Analysis  of  Purnomo  Cl] 
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Figure  10.  Heat  Transfer  Rate  (Q)  vs.  Condenser 
Length 
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APPENDIX  A 


USER'S  MANUAL 

This  appendix  describes  the  data  cards  for  the  use  of 
the  computer  program. 

The  data  is  divided  into  the  COPES/CONMIN  section  and  the 
heat  pipe  analysis  program  section. 

The  COPES  data  is  segmented  into  "blocks''  for  convenience. 
All  formats  are  alphanumeric  for  TITLE  and  END  cards,  F1C  for 
real  data,  and  110  for  integer  data.  Comment  cards  may  be 
inserted  anywhere  in  the  data  deck  prior  to  the  END  card  and 
are  identified  by  a  dollar  sign  ($)  in  column  one.  The  COPES 
data  deck  must  terminate  with  an  end  card  containing  the  word 
'END1  in  columns  1-3. 

Information  is  included  in  this  appendix  pertaining  to 
data  needed  for  single  analysis,  sensitivity  analysis,  or 
optimization  using  COPES/CONMIN. 

UNFORMATTED  DATA  INPUT 

While  the  user's  sheet  defines  COPES  data  in  formatted 
fields  of  ten,  the  data  may  actually  be  read  in  a  simplified 
fashion  by  separating  data  by  commas  or  one  or  more  blanks. 

If  more  than  one  number  is  contained  on  an  unformatted  data 
card,  a  comma  must  appear  somewhere  on  the  card.  If  exponen¬ 
tial  numbers  such  as  2.5+10  are  read  on  an  unformatted  card, 
there  must  be  no  embedded  blanks.  Unformatted  cards  may  be 
intermingled  with  formatted  cards.  Real  numbers  on  an 
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unformatted  card  must  have  a  decimal 
EXAMPLES 


Unformatted  data; 
5,7,1. 3, 1.0  +  20, -5.1 


5,7,1.3,1.0+20, ,-5.1 
5  7  1.3  1.0+20,,  -5.1 
5  7  1.3,  1.0+20  0  -5.1 

Equivalent  formatted  data; 


col-*  10 


20 


30 


40 


50 


DATA  BLOCK  A 


DESCRIPTION:  COPES  TITLE  CARD 

FORMAT:  20A4 

1  2  3  4  S  6  7  8 

Title 


CONTENTS 

Any  80  Character  titJe  may  be  given  on  this  card. 


FIELD 

1-8 
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DATA  BLOCK  B 


DESCRIPTION:  PROGRAM  CONTROL  PARAMETERS 

FORMAT:  7110 

1  2  3  4  5  6  7  8 


NCALC  NDV  NS  V 


FIELD  CONTENTS 

1  NCALC:  Calculation  control 

0-  Read  input  and  stop.  Data  of  blocks  A,B 
and  V  is  required.  Remaining  data  is 
optional . 

1-  One  cycle  through  program.  The  same  as 
executing  ANALIZ  stand-alone.  Data  of 
blocks  A,B  and  V  is  required,  remaining 
data  is  optional. 

2-  Optimization.  Data  of  Blocks  of  A-I  and  V 
is  required.  Remaining  data  is  optional. 

3-  Sensitivity  analysis.  Data  of  blocks, 
A,B,P,Q  and  V  is  required.  Remaining  data 
is  optional. 

2  NDV:  Number  of  independent  design  variables  in 

optimization . 

3  NSV:  Number  of  variables  on  which  sensitivity 

analysis  will  be  performed. 
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DATA  BLOCK  C  Omit  if  NDV=0  in  block  B 
DESCRIPTION:  INTEGER  OPTIMIZATION  CONTROL  PARAMETERS 

FORMAT:  7110 


1  2  3  4  5  6  7  9 


IPRINT  ITMAX  ICNDIR  NSCAL  ITRM  LINOBJ  NACMX1 


FIELD  CONTENTS 

1  IPRINT-  Print  control  used  in  the  optimization 

program  CONMIN. 

0-  No  print  during  optimization. 

1-  Print  initial  and  final  optimization 
information . 

2-  Print  above  plus  objective  function  value 
and  design  variable  values  at  each  iteration. 

3-  Print  above  plus  constraint  values,  direction 
vector  and  move  parameter  at  each  iteration. 

4-  Print  above  plus  gradient  information. 

5-  Print  above  plus  each  proposed  design  vector, 
objective  function  and  constraint  values 
during  the  one  dimensional  search. 

2  ITMAX:  Maximum  number  of  optimization  inter¬ 

actions  allowed.  Default=20. 

3  ICNDIR:  Conjugate  direction  restart  parameter. 

GT.  0-  Scale  design  variables  to  order  of 
magnitude  one  every  NSCAL  iterations.  LT .  0- 
Scale  design  variables  according  to  user  in¬ 
put  scaling  values.  If  not  zero,  NDV  +  1 
is  recommended. 

5  ITRM:  Number  of  consecutive  iterations  which 

must  satisfy  relative  or  absolute  convergence 
criterion  before  optimization  process  is 
terminated.  Default=3. 

LINOBJ:  Linear  objective  function  identifier. 

If  the  optimization  objective  is  known  to  be 
a  linear  function  of  the  design  variables, 
set  LIN0BJ=1.  Def ault=Non-linear . 
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NACMX1:  One  plus  the  maximum  number  of  active 

constraints  anticipated.  Default-NDV=2 . 


f. 
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DATA  BLOCK  D 


Omit  if  NDV=0  in  block  B 


DESCRIPTION:  FLOATING  POINT  OPTIMIZATION  PROGRAM  PARAMETERS 

FORMAT:  7F10 


12345678 


FDCH  FDCHM  CT  CTMIN  CTL  CTLMIN  THETA 


MOTE:  Two  cards  are  read  here. 


FIELD  CONTENTS 

1  FDCH:  Relative  change  in  design  variables  in 

calculating  finite  difference  gradients. 
Default=0 . 01 . 

2  FDCHM:  Minimum  absolute  step  in  finite  dif¬ 

ference  gradient  calculations.  Default= 

0.001. 

3  CT :  Constraint  thickness  parameter.  Default2 

=-0.05. 

4  CTMIN:  Minimum  absolute  value  of  CT  consi¬ 

dered  in  the  optimization  process. 
Default=0 . 004 . 

5  CTL:  Constraint  thickness  parameter  for 

linear  constraints.  Default=-0 . 01 . 

6  CTLMIN:  Minimum  absolute  value  of  CTL  con¬ 

sidered  in  the  optimization  process. 
Default=0 . 001 . 

7  THETA:  Mean  value  of  push-off  factor  in  the 

method  of  feasible  directions.  Default2 

1.0. 
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r 

t 


DATA  BLOCK 

D 

Omit  if  NDV=Q 

in  block  B 

FORMAT : 

4F10 

1 

2 

3  4 

5 

DELFUN  DAB FUN  ALPHAX  AB0BJ1 


I 

l 

FIELD  CONTENTS 

1  DELFUN:  Minimum  relative  change  in  objective 

function  to  indicate  convergence  of  the 
optimization  process.  DEFAULT=0 . 001 . 

2  DABFUN:  Minimum  absolute  change  in  objective 

function  to  indicate  convergence  of  the 
optimization  process.  DEFAULT=0 . 001  times 
the  initial  objective  value. 

3  ALPHAX:  Maximum  fractional  change  in  any  design 

variable  for  first  estimate  of  the  step  in 
the  one-dimensional  search.  DEFAULT =0.1. 

4  AB0BJ1:  Expected  fractional  change  in  the 

objective  function  for  first  estimate  of  the 
step  in  the  one-dimensional  search. 

DE FAULT =0 . 1 . 


REMARKS 

1)  The  DEFAULT  values  for  these  parameters  usually  work  well. 
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DATA  BLOCK  E  Omit  if  NDV=0  in  block  B 

DESCRIPTION:  TOTAL  NUMBER  OF  DESIGN  VARIABLES,  DESIGN 
OBJECTIVE  IDENTIFICATION  AND  SIGN 

FORMAT:  2110, F10 

12345678 


NDVTOT  IOBJ  SGNOPT 


FIELD  CONTENTS 

1  NDVTOT:  Total  number  of  variables  linked  to  the 

design  variables.  This  option  allows  two  or 
more  parameters  to  be  assigned  to  a  single 
design  variable.  The  value  of  each  para¬ 
meter  is  the  value  of  the  design  variable 
times  a  multiplier,  which  may  be  different 
for  each  parameter.  DEFAULT=NDV. 

2  IOBJ:  Global  variable  location  associated  with 

the  objective  function  in  optimization 

3  SGNOPT:  Sign  used  to  identify  whether  function 

is  to  be  maximized  or  minimized.  +1.0  indi¬ 

cates  maximization.  -1.0  indicates  minimi¬ 
zation.  If  SGNOPT  is  not  unity  in  magnitude, 
it  acts  as  a  multiplier  as  well,  to  scale  the 
magnitude  of  the  objective. 
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DATA  BLOCK  F  Omit  if  NDV=0  in  block  B 

DESCRIPTION:  DESIGN  VARIABLE  BOUNDS,  INITIAL  VALUES  AND 
SCALING  FACTORS 

FORMAT:  4F10 

12345678 


VLB  VUB  X  SCAL 


NOTE:  READ  ONE  CARD  FOR  EACH  OF  THE  NDV  INDEPENDENT  DESIGN 

VARIABLES. 

FIELD  CONTENTS 

1  VL3:  Lower  bound  on  the  design  variable.  If 

VLB . LT . -1 .OE+15 ,  no  lower  bound. 

2  VUB:  Upper  bound  on  the  design  variable.  If 

VUB .GT . 10 .E+15  ,  no  upper  bound. 

3  X:  Initial  value  of  the  design  variable.  If 

X  is  non-zero,  this  will  supercede  the  value 
initialized  by  the  user-supplied  subroutine 
ANALIZ. 

4  SCAL:  Design  variable  scale  factor.  Not  used 

if  NSCAL.GE.O  in  BLOCK  C. 
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DATA  BLOCK  G  Omit  if  NDV  =  0  in  block  B 
DESCRIPTION:  DESIGN  VARIABLE  IDENTIFICATION 

FORMAT :  2110, F10 

12345678 


NDSGN  IDSGN  AMULT 


NOTE:  READ  ONE  CARD  FOR  EACH  OF  THE  NDVTOT  DESIGN  VARIABLES 
IN  THE  SAME  ORDER  AS  IN  BLOCK  F. 


FIELD  CONTENTS 

1  NDSGN:  Design  variable  number  associated  with 

this  variable. 

2  IDSGN:  Global  variable  number  associated  with 

this  variable. 

3  AMULT:  Constant  multiplier  on  this  variable. 

The  value  of  the  variable  will  be  the  value 
of  the  design  variable,  NDSGN,  times  AMULT. 
DEFAULT=1 . 0 . 
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DATA  BLOCK  H 


Omit  IF=NDV  0  in  block  B 


DESCRIPTION:  NUMBER  OF  CONSTRAINED  PARAMETERS 
FORMAT:  110 


12345678 


NCONS 


FIELD  CONTENTS 

1  NCONS:  Number  of  constraint  sets  in  the 

optimization  problem. 


REMARKS : 

1)  If  two  or  more  adjacent  parameters  in  the  global  common 
block  have  the  same  limits  imposed,  these  are  part  of  the 
same  constraint  set. 


DATA  BLOCK:  I  Omit  if  NDV  =  0  in  bloc*  B,  or  NC0NS=0  in 

block  H 

DESCRIPTION:  CONSTRAINT  IDENTIFICATION  AND  CONSTRAINT  BOUNDS 
FORMAT:  3110 


12345678 


ICON  JCON  LCON 


NOTE:  READ  TWO  CARDS  FOR  EACH  OF  THE  NCONS  CONSTRAINT  SETS. 


FIELD  CONTENTS 

1  ICON:  First  global  number  corresponding  to  the 

constraint  set. 

2  JCON:  Last  global  number  corresponding  to  the 

constraint  set.  DEFAULT=ICON . 

3  LCON:  Linear  constraint  identifier  for  this 

constraint  set.  LCON=l  indicates  linear 
constraints . 
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DATA  BLOCK:  I 


(CONTINUED) 


FORMAT:  4F10 


12345678 


BL  SCAL1  BU  SCAL2 


FIELD  CONTENTS 


1 

2 

3 

4 


SL:  Lower  bound  on  the  constrained  variables. 

If  BL. LT . -1 . OE  +  15  ,  no  lower  bound. 

SCAL1:  Normalization  factor  on  lower  bound. 

DEFAULT= MAX  of  ABS(BL),  0.1. 

BU:  Upper  bound  on  the  constrained  variables. 

If  BU.GT. 1 . OE+15 ,  no  upper  bound. 

SCAL2:  Normalization  factor  on  upper  bound. 

DEFAULT =MAX  of  ABS(BU) ,  0 . 1 . 


REMARKS : 

1)  The  normalization  factor  should  usually  be  defaulted. 

2)  The  constraint  functions  sent  to  CONMIN  are  of  the  form: 

( BL-VALUE ) / SCAL1  . LE .  0.0  and  (VALUE  -  BU)/SCAL2  . LE . 

0.0. 

3)  Each  constrained  parameter  is  converted  to  two  con- 

constraints  in  CONMIN  unless  A3S(3L)  or  ABS(BU)  exceeds 
1. OE+15,  in  which  case  no  constraint  is  created  for  that 
bound. 
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DATA  BLOCK :  P 


Omit  if  NSV=0  in  block  B 


DESCRIPTION:  SENSITIVITY  OBJECTIVES 
FORMAT:  2110 


12345678 


NSOBJ  IPSENS 


NOTE:  TWO  OR  MORE  CARDS  ARE  READ  HERE. 


FIELD  CONTENTS 

1  NSOBJ:  Number  of  separate  objective  functions 

to  be  calculated  as  function  of  the  sensi¬ 
tivity  variables. 

2  IPSENS:  Print  control.  If  IPSENS. GT.O,  de¬ 

tailed  print  will  be  called  at  each  step  in 
the  sensitivity  analysis.  DEFAULT=No 
print . 


DATA  BLOCK:  P 


(CONTINUED) 


DESCRIPTION: 

FORMAT :  8110 

12345678 


NSN1  NSN2  MSN 3  NSN4 


FIELD  CONTENTS 

1-8  NSNI;  Global  variable  number  associated  with 

the  sensitivity  objective  functions. 

REMARKS : 

1)  More  than  eight  sensitivity  objectives  are  allowed.  Add 
data  cards  as  required  to  contain  data. 


* 


DATA  BLOCK:  £  Omit  if  NSV  =  0  in  block  B 
DESCRIPTION:  SENSITIVITY  VARIABLES 

FORMAT:  2110 

12345678 


ISENS  NSENS 


NOTE:  READ  ONE  SET  OF  DATA  FOR  EACH  OF  THE  NSV  SENSITIVITY 

VARIABLES.  TWO  OR  MORE  CARDS  ARE  READ  FOR  EACH  SET 
OF  DATA. 


FIELD  CONTENTS 

1  ISENS:  Global  variable  number  associated  with 

the  sensitivity  variable. 

2  NSENS:  Number  of  values  of  this  sensitivity 

variable  to  be  read  on  the  next  card. 
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DATA  BLOCK:  £  (CONTINUED) 
DESCRIPTION:  8F10 

FORMAT : 

12  3  4 

SNS1  SNS2  SNS3  SNS4 


CONTENTS 

SENSI :  Values  of  the  sensitivity  variable. 

I=1,NSENS.  1=1  corresponds  to  the  nominal 
value . 

REMARKS : 


FIELD 

1-8 


1)  More  than  eight  values  of  the  sensitivity  variable  are 

allowed.  Add  data  cards  as  required  to  contain  the  data. 


1 


DATA  BLOCK:  V 

DESCRIPTION:  COPES  DATA  'END'  CARD 
FORMAT:  3A1 

12345678 

END 


FIELD  CONTENTS 

1  The  word  'END'  in  columns  1-3. 

REMARKS 

1)  This  card  MUST  appear  at  the  end  of  the  COPES  data. 

2)  This  ends  the  COPES  input  data. 

3)  Data  for  the  user-supplied  routine,  ANALIZ,  follows  this. 
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HEAT  PIPE  ANALYSIS 


Data  for  the  heat  pipe  analysis  follows  the  'END'  card 
in  the  COPES  data  deck.  If  the  general  design  capability  of 
COPES/CONMIN  is  not  needed,  the  heat  pipe  analysis  can 
be  run  by  setting  NCALC  =  1  in  field  number  1  of  data  block 
B;  or  in  a  stand-alone  mode  by  using  the  following  main 
program. 


C  MAIN  PROGRAM  FOR  HEAT  PIPE  ANALYSIS 
C  READ,  EXECUTE,  AND  PRINT  RESULTS 
DO  10  ICALC  =  1,3 
10  CALL  ANALIZ  (ICALC) 

STOP 

END 
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DATA  BLOCK:  AA 


DESCRIPTION:  ELEMENT  CONNECTIVITY 
FORMAT:  315 


12345678 


NEL  NSNP  NBAN 


FIELD  CONTENTS 

1  NEL:  Number  of  elements 

2  NSNP:  Number  of  system  nodal  points 

3  NBAN:  System  band  width 
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DATA  BLOCK:  BB 


DESCRIPTION:  ELEMENT  CONNECTIVITY 


FORMAT:  415 


I  =  1,3;  IEL=1 ,  NEL 


1  2  3  4  5  6  7  3 


IEL  ICOR( IEL ,1 ) 


FIELD  CONTENTS 

1  IEL:  The  element  number 

2  IC0R(IEL,I):  System  nodal  point  corresponing 

to  nodal  point  I  of  element  IEL 


REMARKS : 

1)  Number  all  elements  with  convective  boundaries  first  from 
top  to  bottom,  then  number  the  remaining  elements. 

2)  Number  the  nodal  points  of  each  element  moving  in  a 
counter-clockwise  direction. 

3)  The  elements  with  convective  boundaries  have  nodal  points 
1  and  2  located  on  the  convective  boundary. 
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DATA  BLOCK:  CC 


DESCRIPTION:  CONDENSER  GEOMETRY 

FORMAT:  7G10 . 5 


12345678 


CLI  CANGL  RBASEI  R2I  THICKI  BFIN  TZ 


FIELD  CONTENTS 


1 

2 

3 

4 

5 

6 
7 


CLI:  Condenser  length  (inches) 

CANCL:  Cone  half  angle  (degrees) 

RBASEI:  Inside  radius  of  condenser  base 

(inches) 

R2I:  Intermediate  radius  (inches) 

THICKI:  Condenser  wall  thickness  (inches) 

BFIN:  Height  of  fin  (inches) 

TZ:  Nodal  point  temperature  initial  guess 
(degrees  F) 


REMARKS 

])  Set  TZ  equal  to  TSS-Saturation  temperature  of  the 
working  fluid 
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DATA  BLOCK:  DD 


DESCRIPTION:  FINITE  ELEMENT  GEOMETRY 
FORMAT:  515 


12345673 


NDIV 


NEST 


NEFB 


NBOTI  NBOTF 


FIELD  CONTENTS 

1  NDIV:  Number  of  increments  along  the  length 

of  the  condenser. 

2  NEST:  Number  of  the  element  on  the  right  end 

of  the  trough . 

3  NEFB:  The  element  number  with  convective 

boundary  located  at  the  base  of  the  fin. 

4  NBOTI:  The  element  number  with  convective 

boundary  located  at  the  right  hand  of  the 
bottom  side. 

5  NBOTF:  The  element  number  with  convective 

boundary  located  at  the  left  hand  of  the 
bottom  side. 


( 


I 


DATA  BLOCK:  EE 


DESCRIPTION:  DATA  FOR  RUNNING 

FORMAT:  4F10.2 


1  2  3  4  5  6  7  8 


RPM  TSS  TINF  HINF 


FIELD  CONTENTS 

1  RPM:  Rotation  rate  of  heat  pipe  (RPM) 

2  TSS:  Saturation  temperature  of  the  working 

fluid  (degrees  F) 

3  TINF:  Outside  temperature  (degrees  F) 

4  HINF:  Outside  convective  heat  transfer 

coefficient  ( BTU/HR\Fu  °F) 


DATA  BLOCK :  FF 


DESCRIPTION:  CONVERGENCE  CRITERION 
FORMAT:  G10.9 


12345678 


CRIT 


FIELD  CONTENTS 

1 


CRIT:  Convergence  criterion  on  finite  element 

solution  for  incremental  heat  transfer  rate 


DATA  BLOCK:  GG 


DESCRIPTION:  INTERNAL  FIN  GEOMETRY 
FORMAT :  2G10.5 

12345678 

FANGL  20A 


CONTENTS 

FANGL:  Fin  half  angle  (degrees) 

ZOA:  Ratio  of  trough  width  to  fin  base  width 


FIELD 

1 

2 


DATA  BLOCK :  HH 


DESCRIPTION;  INTERNAL  FIN  GEOMETRY 
FORMAT;  15 

1  2  3  4  5  6  7  8 

IFF 


FIELD  CONTENTS 

1  IFF:  (n-1),  where  n  is  the  number  of  rows  of 

the  upper  triangular  fin  section 

NOTE:  See  Figure  A-l. 
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DATA  BLOCK:  II 


DESCRIPTION:  INTERNAL  FIN  GEOMETRY 
FORMAT :  1615  1  =  1, IFF 

1  2  3  4  5  6  7  8 


KFIN(I)  KFF(I) 


CONTENTS 

KFIN:  The  number  of  system  nodal  points  located 

on  the  symmetric  boundary  of  triangular  fin 
section,  but  does  not  include  the  system 
nodal  points  located  at  the  base  of  the  fin 
and  the  apex. 

KFF:  The  number  of  system  nodal  points  located 
along  the  fin  convective  boundary,  but  does 
not  include  the  system  nodal  points  located 
at  the  base  of  the  fin  and  the  apex. 


NOTE:  See  Figure  A-l. 


FIELD 

1 

2 
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DATA  BLOCK:  JJ 


DESCRIPTION:  INTERNAL  FIN  GEOMETRY 
FORMAT :  2G10.5,4I5 

12345673 


DOBF  DOTH  JTC  JLC  JINT  KT 


FIELD  CONTENTS 

1  DOBF:  Number  of  column  within  fin. 

2  DOTH:  Number  of  column  within  though. 

3  JTC:  The  number  of  the  system  nodal  point 

located  at  the  junction  of  the  symmetry 
boundary  and  the  line  of  intersection 
between  the  fin  and  the  condenser  wall. 

4  JLC:  The  number  of  the  system  nodal  point 

located  at  the  center  of  system  coordinates. 

5  JINT:  The  numerical  difference  between  the  two 

adjacent  system  nodal  points  vertically 
at  the  condenser  section. 

6  KT:  The  number  of  rows  within  the  wall  section. 


NOTE:  See  Figure  A-l. 
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IFF 

=  2 

DOBF  '  *  3.0 

JLC  *  17 

KFIN 

a  2  ;  4 

DOTH  a  1.0 

JINT  =  5 

KFF 

*  3  j  6 

JTC  *  7 

KT  a  2 

Figure  A-l. 

Specification  for  Input 
the  Coordinates  of  the 

Data  to  Determine 
System  Nodal  Points 

fcTB 
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APPENDIX  B 


COMPUTER  PROGRAM  LISTING 


**************************** ********** ******** ********** 
* 

ANALYSIS  CF  ROTATING  HEAT  PIPE  ,  USING  TRIANGULAR 
ELEMENT  MODEL 

COMPILED  BY  MAJOR  I GNA TI US. S. PURNQMO  IN  JUNE  1978 


* 

* 

* 

* 

* 

* 

* 

* 

******************************************************** 


MODI  FI  ED  TO  PERMIT  NUMERICAL  OPTIMIZATION 
USING  COPES/CONMIN 

BY  LCDR  WILLIAM  A.  DAVIS .  JR.  IN  SEPTEMBER  1980 


SUBROUTINE  ANALIZ  CHANGES  THE  ARRAY  DF  DESIGN 
VARIABLES  FROM  SINGLE  TO  DOUBLE  PRECISION  AND 
BACK.  COPES/CONMIN  USES  SINGLE  PRECISION  ONLY; 
DOUBLE  PRECISION  IS  MAINTAINED  IN  SUBROUTINE  FUN 
TO  ALLOW  FOR  POSSIBLE  I LL  CONDITIONING. 


10 

20 

30 


40 


£ 

c 

c 

c 

£ 

c 

£ 

£ 

C 

C 


& 


BROUTINE  ANALIZ  (ICALC) 
MMON  /GLCBCM/  ARRAY (750) 
COMMON  /GLCB1/  BAR AY ( 50) 
REAL *8  BARAY 

IF  (ICALC. GT.l)  GO  TO  20 

KRivM-ij!8oo 

CONTINUE 

CONTINUE 

CO  30  1*1.50 

BARA  Y( I ) »D6LE (ARRAY ( I ) ) 

CALL  FUN  (ICALC) 

DO  40  1*1.50 

ARRAY! I ) *SNGL(BARAY(  I ) ) 

CON?  INUE 

RETURN 


END 
GUIDE 

ALFA 

BFIN 

BVIN 

CALFA 

CANGL 

CBASE 

£IXIT 
St  i 

CPHI 

CRIT 


TO  FORTRAN  VARIABLE  NAMES 


FIN  HALF  ANGLE  (RADIANS) 

HEIGHT  CF  FIN  (INCHES) 

HEIGHT  OF  FIN  (FEET) 

COSINE  OF  ALFA 

CONE  HALF  ANGLE  (DEGREES) 

INSIDE  CIRCUMFERENCE  OF  CONDENSER 

.01  cnaiMfTOft  U  SSKSfiS!58 


IhIj!mal'ConOucti  Vt'TY  OF  CONDEf 

mem \  ttasH  \mu> 

COSINE  OF  PHI 
CONVERGENCE  CRITERION 


(FEET) 
EXIT.  (FE^T 


ATE  fIlM 


(BTU/HR) 


********* 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

£ 

c 

c 

c 

c 

c 

c 


c 

c 

£ 

c 

c 

c 

c 

c 

c 

i 

c 

c 

c 

c 

c 

cc 

c 

c 

c 

c 


OIV 

OMTOT 

DOBF 

DOTH 

FANGL 

H 

HFG 
I  EL 
JLC 

JTC 


KFF 
KF  IN 


KT 

NB  AN 

NBOTF 

NBGTI 

NO  IV 

N6FB 

NEL 

NEST 

NSNP 

PHI 

PI 

RBASE 

RBASEI 

REXIT 

iftf4 

THICK 
THICKI 
TPHI 
ZF  IN 
ZGA 


FLOATING  POINT  VALUE  OF  NOIV 
CONDENSATE  MASS  FLOW  RATE 
NUMBER  OF  COLUMN  WIWTHIN  THE  FIN 
NUMBER  OF  COLUMN  WITHIN  THE  TROUGH 
FIN  HALF  ANGLE  (DEGREES) 

CONVECTIVE  HEAT  TRANSFER  COEFFICIENT  (BTU/HR  FT 2) 
LATENT  HEAT  OF  VAPORIZATION  ( BTU/L8M) 

THE  ELEMENT  NUMBER 

NUMBER  OF  SYSTEM  NODAL  PCI  NT  LOCATED  AT 
THE  CENTER  OF  SYSTEM  COORDINATE 
NUMBER  OF  SYSTEM  NODAL  POINT  LOCATED  AT 
THE  JUNCTION  OF  THE  SYNNETRY  BOUNDARY  AND 
ThE  LINE  OF  INTERSECTION  BETWEEN  THE  FIN 
ANO  THE  CONDENSER  WALL 

NUMBER  OF  SYSTEM  NODAL  POINTS  LOCATED  ALONG 
THE  FIN  CONVECTIVE  BOUNDARY 
NUMBER  OF  SYSTEM  NODAL  POINTS  LOCATED  CN  THE 
SYSMMETRIC  BOUNDARY  OF  TRIANGULAR  FIN  SECTION 
NOT  COUNTING  POINTS  AT  BASE  AND  APEX 
NUMBER  OF  ROWS  WITHIN  THE  WALL  SECTION 
SYSTEM  BAND  W  IDTH 
LAST  ELEMEMT  AT  BOTTOM  SIDE 
FIRST  ELEMENT  AT  BOTTOM  SIDE 
NUMBER  OF  INCREMENT 
ELEMENT  NUMBER  AT  BASE  OF  FIN 
NUMBER  OF  ELEMENTS 
ELEMENT  NUMBER  AT  END  OF  TROUGH 
NUMBER  OF  SYSTEM  NODAL  POINTS 
CONE  HALF  ANGLE  (RADIANS) 

PI 

INSIDE  RADIUS  OF  CONDENSER  BASE  (FEET) 

INSIDE  RADIUS  OF  CONDENSER  BASE  (INCHES) 

INSIDE  RADIUS  OF  CGNCENSER  EXIT  (FEET) 

ifSi  ft 

CONDENSER  WALL  THICKNESS  (FEET) 

CONDENSER  WALL  THICKNESS  (INCHES) 

TANGENT  OF  PHI 
NUMBER  OF  FINS 

RATIO  OF  TROUGH  WIOTH  TO  FIN  BASE  WIDTH 


SUBROUTINE  FUN  READS  INPUT  DAT Af  PERFORMS  HEAT 
TRANSFER  ANALYSIS,  ANO  PRINTS  RESULTS. 


92 


ooo  ooo  non  onnn  ooonono 


SUBROUTINE  FUN  (ICALC) 

IMPLICIT  REAL*8(A-H.0-Z) 

COMMON  /GL0B1/  CL  I , CANGL, R BASE  I , R2 1 , TH I CK I , BF IN, TZ , TSS 
1,TINF,HINF,FANGL,Z0A,ZFIN,B0A,QT0T,RPM 
OIMENS ION  Z( 200),EPS< 200) , HZ (200) , XCOF ( 5 » ,COF ( 5 ) , T { 200 

TEI200'  ’R00W'4' ,Roon  -DELI 200 

COMMON  /APOL/  DOBF ,00TH ,KF INC  50)  ,KFF ( 50)  , I FF , JTC, JLC, J 
1  IN  T ,  KT 

COMMON  /MAF 0/  A ( 200 , 50  I  ,F ( 200 , 1 ) , H (2 00  )  , TS (200  I , T S AT, C 
IK, NEL » NSNP , NBAN, ICQR( 200,3) 

COMMON  /PCRD/  X(200),Y(200),EZERC,8VIV, THICK, TAIFA, APS 
^OENIA^BUZZ)  a(-l.OD0*(Al*ZZ**3/3.0D0+Bl*ZZ**2/2.ODO) 

IF  I CALC*1  READ  INPUT  CATA 
IF  (ICALC. GT.l)  GO  TO  10 


***** 


INPUT  MODE 


***** 


C 

S 


ELEMENT  CONNECTIVITIES 

REAO  ( 5,420)  NEL, NSNP, NBAN 
WRITE  (6,-430)  NEL  ,  NSNP  ,  NBAN 

READ  (5,440)  ( I  EL , (I  COR ( I  EL, l ) , I  =  1 , 3 )  ,  I EL*1 , N EL ) 
WRITE  (6,450) 

WRITE (6, 251) (I  EL, ( I  COR ( IEL, I ) , I* 1, 3 ) , I  EL* 1, NEL) 


ThE  CONDENSER  GEOMETRY 


READ  (  5,460  CLI  tC  ANGL  ,RBA  SEI 

WRITE  (6,470  -  “ 

REAO  (5,480) 

WRITE  (6,490 


ULi  ainol  »koa  ici  ,R2I,THICKI,DriiN,i 

d8fF,*w,M* 

)  NDIV,NEST,NEFB,NBOTI ,NBOTF 


BFIN.TZ 
-  n 
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OATA  FOR  RUNNING 

REAO  (5,500)  RPM,TSS,TINF,HINF 
WRITE  (6,510)  RPM , TSS , TI NF  ,H I NF 

THE  CONVERGENCE  CRITERI AN 

READ  (5,520)  CRIT 
WRITE  (6,530)  CRIT 

INTERNAL  FIN  GEOMETRY 

REAO  (5,540)  FANGL , ZO A 
WRITE  (6,550)  FANGL, ZOA 
READ  (5,560)  IFF 
WRITE  (6,570)  IFF 

REAO  (5,580)  (KFIN(I),KFF(I),I*1.IFF) 

REAO  (5,590)  DOBF, DOTH, JTC , JLC, J INT,KT 
NH8*NEFB/2 

NTM*NBOTI+(NBOTF-NBOTI ) /2 
NBF«NB0TF*1 

WRITE  (6,600)  ICOR ( NBOTI ,2), I  COR (NEF8, 1), IC0R(NTM,2)» I 
1C0R(  NEST ,11,  ICOR(  NBQTF  ,  l) 

RETURN 

IF  I  CALC *2  PERFORM  THE  HEAT  TRANSFER  ANALYSIS 
IF  (ICALC. GT. 2)  GO  TO  360 


***** 


EXECUTION  MODE 


***** 
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CONVERT  UNITS  OF  ALL  DIMENSIONAL  PARAMETERS 

TO  FEET.  CONVERT  UNITS  OF  ANGLES  TO  RADIANS. 

CL-CL1/12.0D0 

R2*R  21/12. COO 

RBAS£*R8ASEI/ 12.000 

eVIN»8FIN/12.0D0 

01  V*0F  L0AT(  NO  I V ) 

PI*3. 1415926535 897900 
PH I *2. ODQ*CANGL*PI /3 60.0 DO 
SPHI*DS  IN  I  PH  I  ) 

CPHI»OCCS(PHI) 

TPHI  *OTAN(  PHI ) 

CELX-CL/CIV 
C§ASE*2.Q00*PI*RBASE 
REXIT=RBASE+CL*TPHI 
CEXIT=2.000*P  I*REX IT 
THICK*THICKI / 12.000 
ALFA*FANGL*2.0D0*PI/360.  0D0 
SALF A=0S I N ( ALFA) 

CALFA=QCQS ( ALFA  I 
TALFA=DT  AN  (ALFA) 

EZER0*2. QC0*8VI N*T  AL  FA 

BOUNDARY  CONDITIONS  AND  TEMPERATURE  ESTIMATES 
ALONG  THE  FIN  BOUNDARY 

DO  20  NT I NF*N80TI » NBGT  F 

TS(NTINF)*TINF 

DO  30  NNT*N8F,NEL 

TS (NNT ) *0.000 

H(  NNTI  *0.  ODO 

DO  40  IGT»1*NEST 

IE*ICOR ( IGT  »2 ) 

T(  IE  )*TZ 
IG*IC0R(NEST, 1) 

Cm4gA»r4m*2 .0D0*P 1*60.000 

DO  50  KL*N8CT I #  NBOT F 
HI KL )=HINF 


FI  FN= 
TSAT: 


=H  INF 
■TSS 


EPSO»ZOA*EZERQ 
E0A*BV I N/ ( EZERO/2 .000 ) 

ZF IN*C  BASE/ ( EZER0+  E  PSC ) 

SURFAR*ZF  IN* (2.000*( BV  IN /CAL FA l+EPSO ) 
EPSEX* ( CEXIT-( ZFIN*EZERO ) J/ZFIN 
BETA*< EPSEX-EPSOi/DI V 
ZZERO BV  IN/ CAL FA 
ZA *0.000 

DO  60  NSAT* 1 tNEST 
TSINSAT )»TSAT 

ts?^o4ISayt«1nn^^-?P?lm 

QT-O.ODO 

QBTOT* 0. CDO 

QT  1*0.000 

OTF*C. 000 

0 TOT-O. ODO 

DMTQT*0 .000 

NK*NOI  V*1 

00  350  NI*I*NK 

RiNIi»R2*NI*DELX*SPHI 

EPS(NI1«EPSC*NI*BETA 

APS* EPS (N  I ) 
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NCCAL  POINT  COORDINATES 

CALL  COORO 
ZU)=ZA 

00  70  IZEL=1,NEFB 

SS:ia8Hiit:H 

XE*X( NA)-X(NB ) 

Y£-Y(NA)-Y(NB) 

£LZ=0SQRT(XE**2+YE**2) 

)  Zl IZEL+1 )*Z ( IZELJ+ELZ 

XZB=X  ( ICCR (NHB* 1 )  )  -X  (  I  COR  (  I»  2  )  ) 

YZB=Y( ICCR(NHB.l) ) -Y ( I  COR ( I » 2 ) ) 

ZB*0SQRT (XZB**2+YZB**2) 

ZC-ZZERO 
I  M=  1 

PARABOLIC  TEMPERATURE  DISTRIBUTION  ALONG  THE  FIN 
BOUNOARY, USING  LAGRANGE  INTERPOLATION 

)  TP1*T (ICCR (1,2) ) 

TP2*T(  ICOR  (NHBf  I)  ) 

TP3*T(ICCfi<N£FB,l ) ) 

AP 1=  TP l/l ZB* ZC ) 

AP2=TP2/ ( ZB*( ZB-ZC ) ) 

AP3*TP3/ ( ZC* ( ZC-ZB1 ) 

BPl*-(  ZB+ZC)*APi 
BP2»-ZC*AP2 
BP  3*— ZB*AP3 

TC*0. COO 
CO  90  NY-1, NEST 
)  TC-TC+T  { ICCR  l  NY  »2 )  ) 

AY*0FL0AT(NY+1) 

TF*(TC+T< ICCR (NY, I) J  +  AY*T  S (NY ) ) / (2.000*AY) 

SOL  10-FLUID  PROPERTIES 
HFG- 1097.200-0. 601 8750 0*TS <11 

RHOFINI 1*62. 77400-0. 0025569800*TF-0. 00005357200*TF**2 
CF(NI) *0.303400+0. 000738927DO*TF-0.00000147321DO*TF**2 
UF (N 1 1 *0 .00 1397D0-0. 0000 1466900*  TF  +  G.  00000006  31253 00* T 
lF**2-0 .0000000000976 56 9D0*TF**3 
LF(NI) -36CC*UF(NI 5 
CM (NI 1-231. 7772D0-0. 02 222DO*TSOL ID 
CM ( NI 1 -8. 776 +0.0026500*7 SO LID 
CM (NI 1*1.0 
CW (N 1 1*20000 .0 


SiSVi 


RHOF ( NI )**2*OMEGA**2*HFG*CPHI*CALFA*R( NI) 


INITIAL  FILM  THICKNESS 
0EL(1)*0. 0000675200 
IF  (NI.GT.l)  GO  TO  100 

OEL( II  *  1. 1C7*( ( ( TSAT- TI NF ) *CF ( NI ) / { UF ( NI ) *HFG*3600 .000 
1  I )  **0 .25)*((UF(NU/(RH0F(NI)*0MEGA)J**0.51 
100  CONTINUE 

AVERAGE  ELEMENT  CONVECTIVE  COEFFICIENT  ALONG 
THE  FIN  BOUNOARY 


ZSTAR* ZZERC-QEL (NI  l/CALFA 
AZZ*OEL( Nil /SALFA 
ZZ*  ZSTAR 

AZS*DABS(4*CF(NI)*UF(NI)*HDEN(A1 ,B1,  ZZ )/ CONST  >**0.2500 
HAC-O. 000 
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OPTIMIZATION  OF  AN  INTERNALLY  FINNED  R0TATIN6  HEAT  PIPE.(U) 
SEP  60  W  A  DAVIS 


non  non  non 


110 

120 


130 

140 


00  190  I EL*1 t NEFB 

az-zued 

BZ*Z(IEL+1) 

IF  (ZSTAR.LE.BZ)  GO  TO  110 
GO  TO  120 

IF  (HAC.NE.0.0D0)  GO  TO  180 
BZ-ZSTAR 

IF  (  IEL.NE.1 )  GO  TO  130 
AK*(8Z-AZ) /5.0D0 
ZZ*  AK 
GO  TO  140 
AK*(BZ-AZ)/4.0D0 
ZZ*AZ 
ZEL*4*AK 
00  1  50  NH*1»  5 
FZ (NH) «CABS ( CF(N I ) **3 *C0NST/ ( 4*UF  ( NI ) *HDEN( A 1 , 
1**0.2500 
150  -ZZ-ZZ+AK 

CGNH*AK* {HZ < 1)+4*HZ( 2) +2*HZ( 3) +4*HZ<  4) +HZ ( 5) )/ 
IF  (ZSTAR.E0.8Z)  GO  TO  160 
H( IEL)*CCNH 
GO  TO  190 
160  AZ=ZSTAR 

HAZ*CCNH*<AZ-Z(  I  EL  )  ) 

0£LA*AZS 
170  BZ*Z ( I EL+ 1) 

DELB*iBZ-ZSTAR)*AZZ/ ( ZZERG-ZST  AR ) 

D£LZ=( OELA+OELB 1/2.000 
FAC* (BZ-AZ1*CF{NI)/DELZ 
H(IEL)*(HAZ+HAC>/<  8Z-ZIIEU) 

GO  TO  190 
180  AZ=Z { I  EL ) 

0ELA=DELB 
FAZ* 0.000 
GO  TO  170 
190  CONTINUE 

NET I*NEFB+l 
DO  200  IEL»NETI,NEST 
200  F( IEL)»CF(NI)/0EL(NI ) 

ENTRY  INTO  THE  FINITE  ELEMENT  SOLUTION 

CALL  FCRMAF 

CALL  BANOEC  (NSNP, NB AN , 1 ) 

THE  TEMPERATURE  DISTRIBUTION 


210 


220 


)  ) 


CO  210  NT* 1, NSNP 
T( NT)*F( NT  fl ) 

TI  B{  NI  i*T( ICOR( NBOTI *  2 
TTTNI)*T(ICCRINEFB, 1) ) 
TBMtNI  )*TUC0R(NTM,2)  ) 
TE(NI)*T(ICOR(NEST, 1 ) ) 
TB<NI)*T( ICCRlNBOTF  *  1 ) ) 
TTS*Q.  OOC 
“  T20  NS« 


CO  21 


* 1» NSNP 


TTS*TT  S*T ( NS ) 

BoKbSHWtt 


Q  AT  THE  BOTTOM  SIDE 
QBI*Q.  OOC 

00  230  IEEL*NB0TI,N80TF 
NKA*ICOR ( I  BEL *1 ) 

HSlftMiffttii 


B  If  ZZ  )  )  ) 
( 3*ZEL) 
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230 


240 


250 


270 


YB*Y(NKA)-Y(NKB) 

ELB»DS QRT(X8**2+YB**2) 

QBI»QBr+(T{NKA7+T{NKBl-2*TS{  IBEL  )  )*ELB*H(  I  BEL  1/2.0  00 
QB(NI)*GBI*0ELX 

ITERATION  UNTIL  CONVERGENCE  CRITERIA  IS  MET 

IF  ( IM.EQ. 1 )  GO  TO  240 
QJ*Q8I 
GO  TO  250 
Q1 *QBI 
I  M=2 

MfccHWo  TO  260 
Ql=QJ 
GO  TO  80 


260  m^c?;sagmfHFG 

Cl*RHCF ( N I i**2*0MEG A**  2*R (NI)*SPHI/( 3*UF ( NI ) ) 
XC0F ( I ) a-DFTCT 
XC0F (2 )®0 .CC0 
XCOF(3)*O.ODO 
XC0F(4I=C1*EPS<NI ) 

XCOFI5  )*C1*T  ALFA 
M*4 

CALL  OPQLRT  (  XCOF , COF  .  M*  ROCTR,  ROOTI ,  I  ER ) 

IF  (RGOTR(l) .GT.0.0D0)  GO  TO  270 
IF  (R00TR(2) .GT. 0.000)  GO  TO  280 
IF  (ROOTR13J .GT.0.0D0)  GO  TO  290 
IF  (RG0TPI4) .GT. 0.000)  GO  TO  300 
WRITE  (6*610) 

WRITE  (6*620)  ( ROOTR  (  I)  ,  I»  1,  4 ) 

GO  TO  6100 


THE  CONDENSATE  THICKNESS 


SSL}SI?fJ*RQ0TR(l, 

280  OEL(NI+l)*RQOTR<2! 


290 

300 

310 


-  'RCCTR13) 
GO  TO  310 
0ELIN241) 

QEL=0. OOC 
IF  (NI.NE.l) 


«R00TR(4) 


320 


330 


GO  TO  320 

Q  FROM  THE  TOP  SIDE 

Q  THROUGH  FIN 

00  330  ICEL*1,3 
KA»ICOR(ICEL,l> 

KB»ICOR( ICEL*  2) 

XQEL*X ( KA) -X  <  KB ) 

YQEL»Y(KBi-Y(KA) 

ELM»DSQRT  (XQEL**2+YQ  EL  **2  ) 

QEL*QEL+12*TS(IQEL)-T(KA)-T(KB) )*ELM*H(IQEL)/2.000 
CONTINUE 

CINC (NI ) *QEL*DELX 
AMT0T(NI )«DMT0T 

SfIg^£ls?5£LX*ZFIN*2 

8tsmfi2:2F,N*2 


G  THROUGH  TROUGH 
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00  340  IQEL»4,4 
KA*IC0R( I  CEL* 1) 

K8*IC0ft< I QEL ,2) 

XQEL»X (KA)-X(KB) 

YQEL*Y (K8) -Y (KA) 

ELM*DSQRT( XQEL**2+YQEL**2) 

GTRF* ( 2*TS (  IQEL )-T (KA )-T( KB) ) *ELM*H(  I3ELJ/2.0D0 
340  CONTINUE 

StAt^AL1? n1T*CINC?nV f  +  QT  INC(NI) 
QTRFT=QTRF*OELX*ZFI  N*2. 

CTFaQTF+fiTRFT 
350  CONTINUE 
RETURN 

360  CONTINUE 


*****  CUTPUT  MODE 
WRITE  16,630) 

CO  370  NR*i,NDIV 
WRITE  (6,640)  NR,  Q I 


***** 


LU  J/U  NH* ItNUlV 

370  WRITE  (6,640)  NR, Q I NC ( NR ) , QT I NC( NR ) , 37 OT AL (NR ) 
wR ITE  (6,650)  QT,  Q  TF 

WRITE  (6,660)  CL  I , CANGL, RBASE I , R 21 , TH I CK I , BF I N ,RPM , TSS 
1,TINF,HINF,CRIT,FANGL,ZQA, IFF 
WRITE  (6,670)  BOA, ZOA, 2FI N,BVIN, SURFAR 
WRITE  (6  ,680 ) 

00  380  NP»1,NSNP 

380  WRITE  (6,690)  NP, X ( NP  )  , Y( NP ) , T(NP ) 

WRITE  (6,700) 

CO  350  KKL=1, NBOTF 
NKX* ICOR (KKL , 1 ) 

NKY*IC0R(KKL,2) 

XP*X(NKX)-X(NKY) 

YP*Y  (NKX  )-Y  ( NKY  l 
EXY*OSORT ( XP**2*YP**2 ) 

CEP* DA BS ((T(NKX)+T(NKY)-2*TS(KKL) ) *E  XY  *H ( KKL ) /2.00  0) 
CEP*GEP*CELX 

390  WRITE  (6,710)  KKL,H( KKL) ,EXY ,QEP 
WRITE  (6,720)  CRIT 

WRITE  (6,730)  HFG , ZF I N,H { NBOTF ) , TS AT , R PM , QTOT , QT,  FANGL 
WRITE  (6,740) 

00  400  NR*1 , NOIV 

400  WRITE  (6,750  NR,  DEL  (  NR)  ,  QB(  NR)  ,  AMTOT  (  NR )  ,T  I B  (  NR)  ,  TT  ( N 
1R),TE(NR),T8(NR) 

WRITE  (6.760) 

00  410  NG=1»NDIV»2 

410  WRITE  (6,770)  NG,CW(NG),CF(NG) ,  RHOF  (  NG  ) ,  UF  ( NG  )  ,  EP  S(  NG  ) 
1,R(NG)  ,TBM(NG) ,QINC(NG) 

RETURN 

420  FORMAT  (315) 

430  FORMAT  ( /2X, 15HN0 .0 F . ELEMENT $*, l 5 , 10X , 34HN0. OF . SYSTEM 
440*  FORMAT  (415.  —  »  •  l«  -  BANDED-, .5. 


a  »EI2.5,/,4X,4HTZ*  ,612.5) 

480  FORMAT  (515) 

490  FORMAT  (4X.6HNDIV*  , I 1 0 ,/ , 4X ,6HNEST*  ,110,/, 
1  ,  I10,/,4X,7HNB0TI*  ,110,/  ,4X  ,7H  NBOTF*  ,U0) 
500  FORMAT  (4F10.2) 

510  FORMAT  (4X,5HRPM*  , E12.5,/ ,4X  ,5HTSS*  ,E12.5, 


4X, 6HN  E  FB! 


f  ru  r\n  m  i 

)  FORMAT  (4X,5HRPM*  , E12.5,/ ,4X  ,5HTSS*  , El 2  .5 , / ,4X, 6HT IN 

i1F0RAS)2(Si^:lf’6HH,NF’  ’E12-S1 
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530  FORMAT  {  4X»  6HCR  IT*  ,E12.  5) 

540  FORMAT  (2G10.5) 

550  FORMAT  i 4X»  7HFANGL*  ,  £12 . 5 ,/ ,4X, 5HZQA*  ,E12.5) 

560  FORMAT  (15) 

570  FORMAT  (4X,5HIFF*  ,110) 

188  n»sr 

600  FORMAT  ( / //5X  ,4HTI B * , I  5 , 10X, 3HTT * , 1 5 , / , 5X ,4HT BM*, 1 5 , 10 
IX, 3HTE=, 15,/ »6X»3HTB* ,15) 

610  FORMAT  (// 10X , 17HCRASH, CRASH,  CRASH) 

620  FORMAT  ( //5X,4( E12. 7 , 3X) ) 

630  FORMAT  < 2X, 7HEL EMENT , 2 X , 4HQF IN, 1 7X , 7HQ TROUGH , 15X, 6HQT0 
1TAL) 

640  FORMAT  ( 4X, I 5,E 12. 5, 10X,E12. 5 , 10X, E12. 5) 

650  FORMAT  <//,4X, 11HQF IN  TOTAL* , E 12 . 5, 10X , 15HQTRCUGH  TOTA 
1L*  ,E12.5) 

660  FORMAT  (/////, 4X*  5HCL I  *  , E  12. 5 »  5 X, 7HCANGL*  ,E12.5,/,4X 
1  ,8HRBASE I  *  ,E12.5,2X,5HR2I*  , E12 .5, /, 4X, 8HTHICKI*  ,E1 2 
2.5,2X,6HBFIN=  ,E12.5,/,4X, 5HRPM*  ,  El  2 . 5 , 5X  ,  5HTSS*  ,E12 
3.5,7,4X,6hTINF=  , E12 . 5, 4X, 6HHINF*  ,E12.5,/,4X  ,6HCRI T* 
4.E12.5 ,4X,7HFANGL*  ,  El 2 .5  , /  ,4X  ,5  HZOA*  , E12.5, 5X, 5HI FF= 
5  ,110) 

670  FORMAT  ( lhl , //2X.4HB0 A*, G 1 2. 5. 5X. 4HZ0 A*, G 12. 5 , 5X, 5HZF I 
1N*,G12.5,5X,5H8FI N=,G12.5,5X ,13HSURFACE  AREA*,G12.5) 
680  FORMAT  ( //5X,2HNP, 6X,  1HX, 12X, 1HY, 12X, IHT) 

690  FORMAT  l /2X  ,  13 , 3X, 3 ( F 10  .6, 3X ) ) 

700  FORMAT  ( //2X ,2HEL , 8X , 1HH , 1 IX ,9HEL-LENGTH , 15X , 4HQ-EL ) 
710  FORMAT  ( / 2X,  12, 3X, E 12 . 5, 3X , E 12 . 5 , 10X , E 12. 5) 

720  FORMAT  < /2X , 22HCCNV ERGENCE  CR ITER  IAN* , E15 .8 ) 

730  FORMAT  ( 1H  , // ,  5X , 4HHFG* , E 12. 5 ,/ , 5X , 1 1HN0.CF  FINS*,E12 
1.5,/ »  5X , 6FH-0UT  =  , E12.5,/«5X, 5HT SAT* , E 1 2. 5, /, 5X, 4HRPM* , 
2E12.5,/,5X,6HQ-80T*,E12.5,/,5X,6HQFIN  * , E12 .5 , / , 5X , 11H 
3FALF— ANGLE*,  F8.  31 

740  FORMAT  ( 1H0 , 6X ,  1HJ,  4X  ,  14HF ILM  THICKNESS, 6X, 8HQ-INCREM, 
16X,8HMASS-TCT,7X,3HTIB,8X,2HTT,10X,2HTE,8X,2HTB) 

750  FORMAT  ( 1H  , 4X, 14,  4X, F 12. 10, 4X,F 10.4, 6X,F 9. 5 , 6X ,F 5. 1 ,6 

IX. F5.1 ,6X,F5.1,6X,F5.1 ) 

760  FORMAT  ( 1HC, 6X, 1H J , 6X, 6HK- WALL ,4X ,6HK-F I LM.3X .7HDENS I T 

IY. 4X »9HV ISC— FILM,6X,7HEPSI LON, 5X , 6HR AD IUS*5X»3HTBM,5X» 
25HQ-BCT) 

770  FORMAT  ( 1H  , 4X, 14, 4X, F7. 3, AX,F 6. 4, 4X , F 6. 3 ,4X , F 9.7 ,4X, F 
19.7,4X,F7.5,5X,F5.1, IX, 1PIC12.3) 

END 

SUBROUTINE  COORD 
IMPLICIT  REAL*8(A-H.C-Z) 

COMMON  /GLOB  1/  CLI ,C ANGL , RBA SE I , R2 1 , TH I CKI , BF I N,T Z , TSS 
1 ,T INF, HINF»FANGL,ZuA,ZFIN, BOA , QT OT , RPM 
COMMON  /PCRO/  XI200)  , Y (200 ) , EZERO, BVIN, THICK, TALFA, APS 
COMMON  /APGL/  DOBF , DOTH, KF I N( 50) ,KFF ( 50) , IF F , JTC, JLC , J 
1INT.KT 

0fcLH*8  VI N/0C8F 
X(  11*0. OCO 
Y(1)*THICK+8VIN 
N*1 

00  20  1*1, IFF 
ICA*KF I N ( I ) 

ICB*KFF( I ) 

CBA-DF LuAT ( ICB-ICA) 

AN *0.000 
00  10  1 1* ICA,  ICB 
X( 1 1 )*X( 1 ) ♦N*AN*DELH*T  ALFA/CBA 
Y( 1 1 )*  Y( 1)-N*DELH 
10  AN*AN+1.0C0 
20  N*N+i 
AN-0.000 


^0°5iCj«jfct JLC,JINT 
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Xt 

Yl J) «( 1. COG-AN/DOTH ) *  THI CK 
DO  30  JJ*1, ICD 

X( J+JJ1«X(J)+JJ*EZER0/ (2*t CBA+l .000 ) ) 

30  Y(  J+ JJ  )*  Y{  J ) 

CO  40  K*1,KT 

X<  J+JJ+Ki*X(  J+JJH>K*APS/(2.0D0*KT) 

40  Y(  J+JJ+K)*YI  J) 

50  AN* AN+l.OCO 
RETURN 
END 

SUBROUTINE  FCRMAF 
IMPLICIT  REAL*8(A-H,0-Z) 

£(JmMQN  1?GlB^  l)  *  ElF  ! i AN&L » R^AS  EI»R2I,THICKI»BFIN»TZ»TSS 
l » TINF,  HINF ,FANGL, ZOA  »ZFIN,BOA,QTCT»kPM 
COMMON  /PCRD/  X (200 > , Y ( 200 >, EZERO, BY  IN . THICK , TALFA , APS 
COMMON  /MAFG/  A( 200 » 50 ) » F <200 , 1) , H(200 I , TS <200 ) ,T S AT , C 
lKtNELfNSNPtNBAN* ICOR (200,3) 

DO  20  N*I » NSNP 
F(N» 1) *0*  CDC 
00  10  MA*1,NBAN 
10  AIN.MA) *0.000 
20  CONTINUE 

CO  70  IEl*L,NEL 
IA*IC0RI  I  EL  ,  I) 

IB*I C0R1 1  EL, 2) 

IC*ICCR(IE(.,3) 

8(  1I*Y(I8  )-Y(  IC> 

6  (  2  )  *  Y (  I C  )  -Y  < IA) 

B(3)*Y(IaJ— Y(IB) 

C(  li*X( IC)-X4 IB) 

C(2)*X (IA)-X(IC) 

itfteSMMIttUm*.*) 

AS«OABS( (B(l)*C(2J-8<2 J*C(l) J/2.0D0) 

HC*H(IEL)/CK 


JJ*ICCi{7l£L,J) 

KK*ICCR(i£c,K) 


iJ)*C(K)  )/(4*AS) 

*  40 

HEL=HC*EL/6.000 
IF  (J.EQ.3)  GO  TO  40 

IF  (K.EQ.31  GO  TO  40 

IF  (J.EQ.K)  GO  TO  30 

EA( J,K)*EA( J,K)+HEL 
GO  TO  40 

30  EA(J,Kl  =*EA  (  J  »K)  *2*H  EL 
40  IF  (KK.LT.JJ)  GO  TO  50 
N'W*KK-JJ*1 

A(  JJ,NW)*A(  JJ,NW)«-EA(  J,K) 

£8St  InuI 

FE*HC*TS< IEL)*EL/2.000 
F< IA , 1 )*F (  IA  » 1 )+FE 
F(IB»1)*F(IB,1I*FE 
CONTINUE 
RETURN 
END 

SUBROUTINE  8AN0EC  (NEQ, MAX8, NVEC ) 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /GLOB  l/  CLI ,CANGL,RBASEI,R21,THICKI ,BFIN,TZ,TSS 
ItT INF*HINF,FANGL, ZOA, ZFIN, BOA.QTOT.RPM 
COMMON  /PCRD/  X ( 200 > , Y (200 ) .EZERO, BY  I N , THICK , T ALF A , APS 
COMMON  /MAFO/  A< 200, 50  ) , F(  200, 1) ,H( 200 ) , T$( 200) ,TSAT,C 


50 

60 


70 


100 


rnn  non  on  o  non  on  nr-non 


1K»  NEL,  NSNPtNBAN* ICOR (200*3) 

L00P«NEQ-1 
DO  20  I-ltLCOP 
MB-I+1 

N8-MIN0I  I+MAXB-1* N  EQ  ) 

00  20  J-HB,N8 
L*  J+2-MB 
D*A(  I,  L)/A  ( I  ,  1 ) 

00  10  MM-1 »NVEC 
)  F(J.MMJ*F(J,MM)-D*F( I, MM) 

MM-MI N0( MAXB-L+l*  NEQ— J  +1 ) 

00  20  K- 1 »MM 
NN-L+K-l 

>  oAJ5o'!il^:0*4<I’NNI 

)  F1NEC«I)*FINEQ»I)/A(NEQ» 1) 

00  50  1-2.,  NEQ 
J-NEQ-IM 

k-minoinec-j+umaxb) 

00  50  MM-1,NVEC 
00  40  L*2,K 
MB-J+L-1 

3  F(  J,MMi*F( J,MM|-A<  J,  L) *F(MB,MM) 

)  F(J»MM)-F<J,MM)/A(J, 1) 

RETURN 

END 

subroutine  cpglrt  computes  the  roots  of  a  real 

POLYNOMIAL  USING  A  NE  Vi  TON-  RA  PH  SC  N  ITERATIVE 
TECHNIQUE. 

SUBROUTINE  OPCLRT  <  XCOF,  COF,  M,  ROOTR,  ROOTI ,  IER  ) 
IMPLICIT  REAL*8  (A-HI,  REAL*8  (O-Z) 

DIMENSION  XCOF(l),  COFi  I),ROOTR{1)»ROOTII  1) 


IF  IT-0 
N-M 
IER*  0 

IF  (XC0F(N*1))  10,40,10 
10  IF  (N)  20,20,60 

SET  ERROR  CCOE  TO  l 


20  IER-1 
30  RETURN 


SET  ERROR  COOE  TO  4 


40 


TO  30 

SET  ERROR  COOE  TO  2 


50  IER- 2 
GO  TO  30 

60  IF  (N-36J  70,70,50 

70  NX-N 
NXX-N+1 
N2- 1 
KJl-N-H 
00  80  L»l»KJi 
MT-KJ1-L*! 

80  COFI MT I-XCCF (L) 

SET  INITIAL  VALUES 


nno  non  non  non  non  on  n 


90  XO*. 00500101 
Y0»0. 01000101 


ZERO  INITIAL  VALUE  COUNTER 

IN*0 
100  X*X0 

INCREMENT  INITIAL  VALUES  AND  COUNTER 

XO*-  10. 0*YO 
YQ— 10.0*X 

SET  X  ANO  Y  TO  CURRENT  VALUE 

X*X0 
Y*Y0 
IN-IN+1 
GO  TO  120 
110  IF  IT* 1 
XPR*X 
YPR*  Y 

EVALUATE  POLYNOMIAL  ANO  DERIVATIVES 

120  ICT*0 
130  UX*0.0 
UY*0 .0 
V-0.0 
YT*Q.O 
XT*1 • 0 
U*C0F(  N+l  J 
IF  (UJ  140,270,140 
140  DO  150  I *1 ,N 
L*N— I* 1 
XT2*X*XT-Y*YT 
YT2*X*YT+Y*XT 
U*U*CCF (L ) *XT  2 
V*V-*-COF  i  L » *YT2 
FI*I 

UX*UX+F  I*XT*C0F(L ) 

UY*UY— FI*YT*CCF ( LI 
XT*XT2 
150  YT*YT2 

SUMSQ«UX*LX+UY*UY 
IF  ( SUMSQ)  160,230,  160 
160  DX*( V*UY-U*UX) /SUMSQ 
X*  X+DX 

0Y»-(U*UY*V*UX ) /SUM SO 
Y*Y*QY 

IF  (  DA8S(  QY ) +DABS( DXI-1.0E-05)  210,170,170 

STEP  I 1ERATI CN  COUNTER 

170  ICT*ICT+1 

IF  ( ICT-5CCI  130,180,180 
180  IF  4  IF  IT  I  210,190,210 
190  IF  (IN-5)  ICO  ,200,200 


SET  ERPCR  CODE  TO  3 


200 

210 


igRx2  30 
00  220  L* 1 ,NXX 
MT-KJ1-L+1 
TEMP*XC0F ( MT) 
XCOF(MT)*CCF( L) 
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280*260,260 


COF{L)»TEFP 

ITEMP-N 

N»NX 

NX*I  TEMP 

IF  (  IF  IT  J  250,  110,250 

IF  ( IF  IT)  240,100,240 

X-XPR 

Y-YPR 

IF  I T*0 

IF  IDABS( Y/XJ-1.0E-04)  280,260,260 

AlPHA*X+X 

SUMSQ»X*X+Y*Y 

N*N-2 

GO  TO  290 

X-0.0 

NX*NX— 1 

NXX*NXX-1 

Y*O.C 

SUMSQ»0.0 

ALPHA»X 

N*N-1 

COF 121 »CCF (2 ) ♦ALPHA*CQ  F( 1 ) 

00  300  L*2  »N 

COFa«-l)=COFlL+l)+ALPHA*C.QF(  L)-SUMS3*C0FIL-1> 

ROOT I ( N2 ) *Y 

RQ0TR(N2)»X 

N2»N2+i 

IF  { SUMSQ)  320,330,320 

Mrm 

IF  (N)  30,30,90 
ENO 
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